System and method for automatic color segmentation and minimum significant response for measurement of fractional localized intensity of cellular compartments

ABSTRACT

A system, a method, and a programmed device for measurement of translocational activity among cellular compartments process magnified images of cellular material exposed to an agent by segmenting and compartmentalizing the images and then measuring fractional localized intensity of two or more compartments in the segmented, compartmentalized image. The measured fractional localized intensities are compared to determine translocation of cellular material among the measured components caused by the agent.

RELATED APPLICATIONS

This application is a continuation of U.S. patent application Ser. No.10/507,442, §371(c)(1),(2),(4) date Mar. 15, 2006, which is a U.S.national phase under §371 of PCT Application PCT/US03/07968, whichclaims priority from U.S. Provisional Patent Application 60/363,889,filed Mar. 13, 2002.

TECHNICAL FIELD

This application concerns assay of biological material by means of highspeed, high throughput cytometry using fractionalized localizedintensities of subcellular components in magnified images.

BACKGROUND ART

Drug discovery screening has historically used simple well-basedread-outs in order to handle a high throughput of lead compounds.However, a given assay currently provides only the information that adrug affects some of the cellular processes that result in the responsemeasured; the exact nature of the target for the drug is not indicated.A cell-based assay is a model system designed to identify compounds thatinteract with a selected molecular target in a specific manner.Cell-based assays are robust in that they approximate physiologicalconditions, and they can yield highly complex information. This requiressophisticated image analysis tools and streamlined data handling.Multi-parameter cell assays, where a response is measured by multiplexedreporter molecules, as well as morphological criteria, have been limitedby the labor-intensive nature of imaging and analyzing subcellulardetails. The power of obtaining more complex information earlier in thescreening process demands effective solutions to this bottleneck.

Dissecting the steps in cellular pathways is important, because multiplepathways converge and diverge to provide redundancy (backup in case ofcellular dysregulation) and as a coordinated response. Whereas a givendrug may result in, e.g., the secretion of a cytokine (such asInterleukin 2 (IL-2) measured as a single parameter response in a wellplate) the researcher does not know which signaling pathway wasutilized, or what other cellular responses were initiated. If thesignaling pathway used also led to cell death, the efficacy of thecandidate drug would be compromised and would fail in costly andcontroversial animal testing. Multiplexed cellular responses need to beinvestigated to eliminate this kind of false positive lead in drugdiscovery.

The complexity of real cell responses also leads to heterogeneitybetween cells, even in a cloned cell line, depending on other factorssuch as cell cycle progression. Thus, a lead which acts on a cellularprocess which is part of DNA synthesis would elicit a response only inthose cells which were in S phase at the time the drug was added. In theclinical situation, continuous infusion can ensure all cells aretreated, and so the lead is a viable candidate. If an average responsefrom all the cells in a well is measured, it may fall below thethreshold for detection and result in a false negative: an effectivedrug is overlooked.

Pharmaceutical companies continually demand faster analysis of screeningtests. Automation has addressed the need to increase data acquisition,but there remain stringent requirements for accuracy, specificity andsensitivity. Preliminary data indicates that the higher the informationcontent of the assay read-out, the less variability there is betweenindividual cells in a responding population. Thus, the total number ofcells needed to attain the confidence level required by the experimentis decreased, resulting in increased throughput. More accurate analysisresults in better dose response information. Higher quality data resultsin better data mining and identification of drugs with significantclinical impact.

Automated quantitative analysis of multiplexed fluorescent reporters ina population of intact cells at subcellular resolution is known.Accurate fluorescence quantification is made possible by technologicaladvances owned by Q3DM, the assignee of this application, and therapidprocessing of this complex image data in turn depends on uniquecomputational processes developed by Q3DM. High throughput microscopyhas only recently become possible with new technological advances inautofocus, lamp stabilization, image segmentation and data management(e.g., U.S. Pat. Nos. 5,548,661, 5,790,692, 5,790,710, 5,856,665,5,932,872, 5,995,143, and U.S. patent application Ser. Nos. 09/703,455and 09/766,390). Accurate, high-speed autofocus has enabled fullyautomated “walk-away” scanning of arbitrarily large scan areas.Quantification is dramatically improved by controlling the intensity ofillumination, and is dependent on accurate autofocus. Advances in imagesegmentation make detection of both dimly and brightly stained cellssimple, so that heterogeneous cell populations can be assayed withstatistically meaningful results. Finally, rapid streaming of highinformation content data depends on efficient image data format andcaching. Together, these technological advances enable retrieval of highquality, image-based experimental results from multiplexed cell basedassays.

SUMMARY AND ADVANTAGES OF THE INVENTION

The most comprehensive tool set that biology drug discovery couldpossess at this time would be one that could integrate technologiesinvolving the acquisition of molecular, cell-structural, andcell-activity data. Application of such a tool set to Biology programswould transform current end-point cell-assays from indirect read-outs(such as IL2) to sub-cellular structural co-localization analyses, andtime-lapse imaging of cell function activities. Quantitative cellimaging technologies will lead the way to the next level of rationaldrug screening and design.

With the automation of this rate-limiting step in the discovery pipelinecomes the possibility of designing high throughput cell-based assaysthat occur in the context of the whole cell, as well as the ability tomultiplex reporter molecules so that complex cellular interactions maybe taken into account. These assays would have measurable economicadvantages, because in obtaining more data earlier, they offer thepossibility of eliminating false positives and false negatives due toheterogeneity in cellular response. The benefit to drug discovery isthat lead compounds are better qualified before they enter costly andcontroversial animal trials.

These improvements and advantages are provided by this invention whichis realized as an assay system and method for automating colorsegmentation and minimum significant response in measurement offractional localized intensity of cellular compartments.

Assay development may yield unexpected results, therefore the inventionprovides the ability to visually inspect cells identified on amulti-parametric data plot. Gating the cells of interest generates animmediate montage on a display screen, and by relocation of the scanningstage, the user can look through a microscope eyepiece for validation.

New assay development demands that the assay system be adaptable tonovel reporter molecules, to new cell types, and to the visual nature ofthe cellular response. To this end, the assay system of the inventionuses object-oriented software framework to be open and extensible sothat features such as hardware, measurement or database functions, andcell classification schemes can be added without having to uncover largeportions of existing software. In particular, adaptive objectrecognition algorithms incorporated into the assay system of theinvention, including image segmentation and tessellation, allow a userto interactively define what constitutes an object of interest in anassay (i.e. the nucleus of a given cell type). Implementedinteractively, these algorithms are incorporated into scanningparameters for that assay, so that automated image segmentation andtessellation can enable the speeds necessary for high throughputscreening.

The accuracy of fluorescence measurements is enhanced due to highperformance features, notably autofocus, lamp stabilization, imagesegmentation, image tessellation, and image measurement algorithms. Theresulting increase in the signal-to-noise ratio of the data is evidentin diminished residuals in a best-fit curve, providing enhancedsensitivity.

BRIEF DESCRIPTION OF DRAWINGS

FIG. 1 is a block diagram of an automated microscopy system according tothis invention.

FIG. 2 is a plot showing a surface of best foci in a 9×15 squaremillimeter area.

FIG. 3 is a diagram demonstrating a process using incremental scanningto obtain a best focus which may be used in the practice of thisinvention.

FIGS. 4 a, 4 b, and 4 c illustrate autofocus circuit dynamic range andsensitivity.

FIG. 5 is a block diagram of an autofocus measurement circuit which maybe used in the practice of this invention.

FIGS. 6 a-6 d, illustrate enhancement of contrast in magnified images bynonlinearly trained (perceptron criterion) 2D image filters.

FIG. 7 is a magnified image of biological material showing tessellationof cell nuclei.

FIG. 8 is a plot showing results obtained using optical feedbackstabilization of a 100 W Hg vapor arc lamp.

FIG. 9 is an illustration of an image table data structure which may beused in the practice of this invention.

FIG. 10 is an illustration of an image table caching algorithm which maybe used in the practice of this invention.

FIG. 11 is a series of three panels of magnified cell images showing theeffect of multichannel image acquisition and processing softwarecomponents on a BIPI-prepared wellplate.

FIG. 12 includes four panels showing how image segmentation andtessellation may be used in the invention to identify and localize cellnuclei.

FIG. 13 illustrates the population statistics for TNFα-induced NFκBactivation in HUVEC cells.

FIG. 14 is a family of curves representing the theoretical dose responseresolution for σ={0.02,0.04,0.06,0.08,0.1,0.12,0.14} (ascending, α=0.05(Type I error probability), β=0.05 (Type II error probability)).

FIG. 15 shows a theoretical dose response curve with μ_(max)=1.0,μ_(min)=0.0, and IC₅₀=8.895×10⁻⁵.

FIG. 16 is a family of plots showing a worst-case response for well-wellmean FLIN response.

FIG. 17 is a family of plots showing a best-case response for well-wellmean FLIN response

FIG. 18 is a family of plots showing the results of a Monte Carlosimulation of dose response regression.

FIGS. 19 a and 19 b are curves that show Monte Carlo confidence intervalwidth versus sample concentrations.

FIG. 20 is a plot showing a translocation response of the G₁subpopulation according to the invention.

FIGS. 21 a and 21 b illustrate inhibition of NFκB nuclear translocationby BIPI compound A.

FIGS. 22 a and 22 b illustrate inhibition of NFκB nuclear translocationby BIPI compound B.

FIGS. 23 a and 23 b illustrate inhibition of NFκB nuclear translocationby BIPI compound C.

FIG. 24 illustrates inhibition of NFκB nuclear translocation by BIPIcompound B in a defined G₁ cell subpopulation

FIG. 25 is a family of curves produced by plotting time required toachieve a minimum significant response as a percent of total NFκBtranslocation on maximum stimulation per well.

FIGS. 26 a and 26 b are illustrations of first and second graphical userinterface (GUI) tools used to initialize a drug screen according to theinvention.

FIG. 27 is an illustration of a third GUI tool used to initialize a drugscreen according to this invention.

FIG. 28 is a flow diagram representing steps executed in initializing adrug screen.

FIG. 29 is a flow diagram representing a drug screen method performed bya programmed computer according to this invention and also representinga software program for causing the method to be performed.

FIGS. 30 a-30 h are geometric figures including nodes, edges, circles,and polygons that illustrate various aspects of tessellation as used inthis invention.

DETAILED DESCRIPTION Fractional Localized Intensity of CellularCompartments (FLIC)

Development of multicompartment models for cellular translocationevents: Many potentially important molecular targets are regulated notonly in their expression levels, but also by their subcellular orspatial localization. In the post-genomics era, illuminating thefunction of genes is rapidly generating new data and overthrowing olddogma. The prevailing picture of the cell is no longer a suspension ofproteins, lipids and ions floating around inside a membrane bag, butinvolves protein complexes attached to architectural scaffolds orchromatin provided by the cytoskeleton, endoplasmic reticulum, Golgiapparatus, ion channels and membrane pores. Cell surface receptors areoriented within the plasma membrane such that they can bind anextracellular molecule and initiate a response in the cytoplasm. Proteincomplexes in the cytoplasm can be dissociated following regulatedproteolysis to release specific DNA binding proteins. These proteinsthen pass through nuclear pores, interact with the chromatinorganization, and regulate gene transcription. Proteins are thentrafficked through the Golgi apparatus where they are readied forfunctionality. Any of these processes can be targets for improvingclinical efficacy and reducing side effects, and as such are importantto understand.

A typical assay to screen for an agonist or antagonist for this processwould measure the movement of a known DNA binding protein from thecomplex into the nucleus. However, by multiplexing reporter molecules,the same assay can also provide information on the receptorinternalization. Thus the user can be sure of which receptor isresponding to the drug when the downstream effect is mediated. Broadlyspeaking, there is a need to visualize protein movement as a response tothe activation of target signaling molecules, specifically receptors onthe external membrane binding their ligand (or drug) and internalizing,as well as transcription factors moving from cytoplasm to nucleus. Thedefinition of discrete compartments has been achieved by compartmentspecific labeling (dyes). The plasma membrane is labeled with alipophilic carbocyanine dye, the cytoplasm with a dye that permeateslive cells but is cleaved in the cytoplasm so that it cannot diffuse outand the nucleus by a membrane-permeant DNA intercalating agent.

Round cells such as lymphocytes present a challenge in identifyingnuclear, cytoplasmic, and membrane compartments. Resolving the sparsecytoplasmic compartment requires high resolution imaging, achieved usinga high numerical aperture objective lens, but this also narrows thedepth of field. There are many subcellular components, organelles orpatterns that also require high-resolution images. Robust autofocusbecomes necessary; this must be done rapidly to meet the demands of highthroughput screening, and to ensure accurate quantification offluorescence intensity. The invention employs an autofocus mechanismthat is an order of magnitude faster than other systems availablebecause it uses a dedicated circuit to measure image sharpness directlyfrom the video signal. Such an autofocus mechanism is robust because itmeasures the change in the optical transfer function (OTF) in a rangethat avoids the contrast reversals inherent in cell images.

The assays used for generating preliminary data will may include asecond reporter molecule that detects a separate protein. Desirably,more parameters may be measured in the same cell over a given timeperiod. By the use of an object-oriented software framework, theplatform employed by this invention is extensible and algorithm-drivenfor maximum flexibility. To give a hypothetical case, a simplecomparison of two spatial patterns (e.g., nuclear and cytoplasmicfluorescence) may be insufficient to determine the efficacy of a lead ona target. Even though the positive response occurs, if this isaccompanied by changes in cell morphology or metabolism, the efficacymay be questionable, or toxicity may be a more important limitation. Theaddition of reporters for these cell properties will require additionalimage analysis algorithms that will classify a set of cellularresponses. Preferably the invention utilizes software pluginarchitecture in order to meet these needs by supporting development ofnew, assay specific processing modules into the image analysiscapabilities of the base platform and existing assays.

A high throughput, microscopy platform utilized in the practice of theinvention is illustrated in FIG. 1. This platform includes a highperformance, research-grade instrument built as an attachment systemaround a commercial fluorescence microscope. The platform includes aprogrammable digital (host) computer 110, preferably a high-end,off-the-shelf Pentium® workstation running the Windows® operatingsystem. Developed and validated assays may be run on this system using astraightforward, point and click operation based upon familiar Windows®motifs. Other components of this platform include a live monitor 115, aninverted Nikon TE300 fluorescence microscope 117 with attachments thatinclude a robotic stage 118, a CCD (charge coupled device) camera 120,and an LED phase contrast light source 125. A stabilized mercury arclamp mechanism 131 provides illumination for operation of the microscopeby way of fiber optic housing 132 and an optical fiber 133. A controlunit 137 drives the stage 118 while it supports a structure 139 (such asa microtiter plate or another equivalent assay device) disposed on thestage 118 in response to commands issued by the host computer 110 andhost computer. An autofocus circuit 142 responds to images of obtainedby the CCD camera 120, providing processed information to the hostcomputer 110. A piezofocus mechanism 148 controls the positioning of thelens mechanism of the microscope 117 in order obtain a best focus. Thehost computer 110 is programmed to perform routines for focusing themicroscope 117, controlling the stabilized mercury arc lamp mechanism131, scanning the structure 139 by moving the stage 118, intaking andstoring sequences of magnified images of biological material on thestructure 139 produced by the CCD 120, processing the magnified imagesby segmentation and tessellation, and performing the translocationprocessing described in more detail below.

Autofocus is critical for scanning large areas due to the variations inslide and coverslip surfaces and mechanical focus instability of themicroscope (particularly thermal expansion) [M Bravo-Zanoguera and J HPrice. Analog autofocus circuit design for scanning microscopy. InProceedings, International Society for Optical Engineering (SPIE),Optical Diagnostics of Biological Fluids and Advanced Techniques inAnalytical Cytology, volume 2982, pages 468-475, 1997]. These effectscombine to effectively create an uneven surface over which the cellsmust be analyzed. An example of this uneven surface is plotted in FIG.2, a 3D plot of a surface of best foci (9×15 mm² area) demonstrating theneed for autofocus. The plot represents the surface of best foci withthe average plane subtracted and shows a range of 6 μm. [MBravo-Zanoguera, B von Massenbach, A Kellner, and J H Price. Highperformance autofocus circuit for biological microscopy. Review ofScientific Instruments, 69(11):3966-3977, 1998]. As can be appreciatedwith reference to this figure, the foci can easily vary through a rangeof ±10 μm over an entire microscope slide. This can have a dramaticeffect on the ability to locate cells. Using a Nikon Fluor 0.75 NA 20×objective, for example, the mean fluorescence intensity has been shownto drop by about 35% at ±10 μm from focus [J H Price and D A Gough.Comparison of digital autofocus functions for phase-contrast andfluorescence scanning microscopy. Cytometry, 16(4):283-297, August1994].

FIG. 3 is diagram demonstrating incremental scanning. Best focus fromthe previous field provides the center of the autofocus search range forthe next field. The degree of focus (sharpness) is then measured forseveral test planes (typically 5-9) at the video interlaced field rateof 60 Hz, and the piezoelectric positioner is moved to best focus. Stagemovement and autofocus require a total of 0.3 s.

Incremental scanning is carried out by moving the stage to a field,stopping the stage, performing autofocus, acquiring the image andrepeating on the next field. This sequence is shown in FIG. 3. If focuslies too close to one extreme of the search range, the center of therange can be repositioned and autofocus repeated.

Autofocus on single fields has been extensively explored and reviewed[Price & Gough]. Further work has extended the techniques to largenumbers of microscope fields and performing high-speed autofocus withreal time image processing [Price & Gough]. Careful attention tomagnification and sampling led to about 100 nm precision (as measured bythe standard deviation) in scans of thousands of microscope fields withfocus achieved in 0.25 s. Further development of this technology led todesign and implementation of an autofocus circuit that tremendouslyreduced the cost of autofocus and improved focus precision to an averageof 56 nm [Bravo-Zanoguera et al.]. Additional theoretical andexperimental studies on the through-focus OTF helped further support thechoice of the correct filter for use in isolating the high frequenciesfor focus measurement [M A Oliva, M Bravo-Zanoguera, and J H Price.Filtering out contrast reversals for microscopy autofocus. AppliedOptics, 38(4):638-646, February 1999].

TABLE 1 Autofocus performance in scanning Results Mean ExperimentalConditions Best Focus Range Combined σ Digital- Area Test Range Max -min (μm) Digital Analog Analog Experiment (mm × mm) Fields Z(μm)Position^(a) Zoom Raw Nonplanar^(b) (μm) (μm) (μm) 1^(c, d) 3.25 × 1.911200 1.465 11 2.0  6.2 2.8 0.146 0.106 0.144 2^(d, c) 3.48 × 2.05 17282.197 11  2.25  7.3 2.4 0.144 0.052 0.088 3^(d, c) 2.60 × 2.05 12962.196 11 1.1  7.9 1.5 0.036 0.027 0.043 4^(d, c) 3.63 × 2.84 2500 2.19611 1.1 16.9 4.1 0.075 0.067 0.006 5^(c, g) 3.87 × 3.04 1296 2.196 11 1.510.5 2.1 0.031 0.027 0.074 6^(h, g) 4.00 × 3.12 1369 2.196 11 1.5 12.81.6 0.028 0.023 0.068 7^(h, i) 3.87 × 3.04 1296 2.196  9 1.5 11.3 2.90.051 0.042 0.122 ^(a)Focus time is (positions + 2)/60 s, or 0.18 s and0.22 s for 9 and 11 test positions respectively ^(b)Plane fit to data bylinear regression and subtracted ^(c)Linear spacing between focus planes^(d)Digital Tracking ^(e)Nonlinear spacing between focus planes of 17,10, 7, 6, 6, 6, 6, 7, 10 and 17 digital units ^(f)48% overlap betweencontiguous fields ^(g)Analog tracking ^(h)Tracked by average of analogand digital ^(i)Nonlinear spacing between focus planes of 22, 10, 7, 6,6, 7, 10 and 22 digital units

FIGS. 4 a, 4 b, and 4 c illustrate autofocus circuit dynamic range andsensitivity demonstrated using (left) a microscope field with no cellsand only a small amount of debris. (Center) shows the plot with autogaindisabled (gain of 1.0), and (right) shows the plot with an autogain of100.

The autofocus circuit illustrated in FIG. 5 was tested under typicalscanning conditions [M Bravo-Zanoguera, B von Massenbach, A Kellner, andJ H Price. High performance autofocus circuit for biological microscopy.Review of Scientific Instruments, 69(11):3966-3977, 1998.]. The purposeof the scanning experiments was to determine the ability of the systemto track focus, as well as measure precision and compare the analogcircuit to the digital system. Seven rectangular areas, ranging from1200 to 2500 fields, were scanned in a raster pattern. The results areshown in Table 1. For the first six experiments, focus was computed fromthe power-weighted average of the focus measurements from 11 differentvertical positions. In the final experiment, this number was reduced to9 vertical positions. Autofocus was performed 20 times at eachmicroscope field. The combined standard deviation (σ) of all autofocustrials (20 fields) in each experiment is shown in Table 1. For everyexperiment, the precision (as measured by σ) was better for the analogcircuit than the digital system. The combined σ s for all sevenexperiments were 0.056 μm (analog) and 0.087 μm (digital). This isalmost an order of magnitude better than the 0.528 μm depth of field(according to the Francon criterion [M Francon, editor. Progress inMicroscopy. Row and Peterson, Evanston, Ill., 1961]) and approaches theminimum vertical step of 0.024 μm (PIFOC piezoelectric objectivepositioner, Polytec PI, Irvine, Calif.).

Table 1 also includes two columns showing the range of best foci foreach area scanned. In one column the raw data range is presented, whilein the next column the range with the average plane subtracted is shown(nonplanar). The raw data range is 6.2-16.9 μm and is largest for thebiggest area. The nonplanar range is 1.6-4.1 μm and is still much largerthan the depth of field of high NA objectives. Other experiments (datanot shown) indicated that the absolute variation from flat increaseswith area as expected. For example, data from larger 10×15 mm² scansrevealed a 6 μm range, and further experience scanning larger areas hasled us to expect as much as a ±10 μm range of foci over an entire slide.This range is even larger (as much as hundreds of microns) for thewellplate formats that dominate many industrial microscopy applications.

The autofocus circuit of FIG. 5 provides very high sensitivity. However,it became clear early in the design process that the analog sensitivitywas much greater than available with subsequent 12 bit A/D conversion.Thus, an autogain circuit was added to augment the dynamic range. FIGS.4 a, 4 b, and 4 c demonstrate the advantage of autogain, and the dynamicrange and sensitivity of the circuit on a microscope field containingonly small debris. We have found that this circuit is so sensitive thatit can track focus even when no cells are present as long as there issomething, even if only minimal debris, in the field.

It would be convenient if simple intensity thresholding were be adequatefor quantitative cell-based assay of cells stained with relativelybright fluorescent dyes. Unfortunately, fluorescently stained cellsinvariably exhibit wide variations in stain because the cell size andfluorochrome densities vary. Accordingly, real-time image segmentationis utilized in this invention for fluorescently stained cellcompartments in order to make subcellular compartment identification andlocalization much less dependent on fluorescence intensity orvariability [J H Price, E A Hunter, and D A Gough. Accuracy of leastsquares designed spatial FIR filters for segmentation of images offluorescence stained cell nuclei. Cytometry, 25(4)303-316, 1996]. Thiswork is illustrated in FIGS. 6 a-6 d, wherein enhancement by nonlinearlytrained (perceptron criterion) 2D image filters to optimally approximatehuman texture classification capabilities is illustrated. In FIG. 6 a,original images of three DAPI stained 3T3 cell nuclei are shown. FIG. 6b illustrates the result obtained with application of conventional imagefilters which do not provide adequate contrast enhancement. In FIG. 6 c,linearly designed 2D filters also do not provide adequate contrastenhancement. FIG. 6 d illustrates that results obtained with nonlinearlydesigned filters provide dramatic contrast enhancement enabling simpleintensity thresholding for a final segmentation step. This methodutilizes nonlinear least-squares optimized image filters to createmarked object-background contrast for automatic histogram-basedthresholding. This contrast makes the final image segmentation much morethreshold-independent, and allows dim cells to be segmented with thesame accuracy as bright ones over a much larger intensity range than waspreviously possible. According to this method, the error function forfilter design is allowed to sum only the error from a predeterminedminimum contrast. Allowing intensities to extend beyond this range,rather than requiring an exact match, substantially improved performance(see FIGS. 6 a-6 d). This method has been tested extensively over manyyears. For example, Price, Hunter and Gough tested ten montage imagescontaining a combined 1,070 fluorescently stained cell nuclei. Eachimage was manually segmented to provide a standard. Then 3×3 through25×25 element 2D filters were trained using the perceptron criterion andnonlinear least squares to optimally approximate the human standard.Each filter was tested against the other 9 images in order to ensurethat the image did not bias this supervised design technique. Very highsegmentation accuracy was achieved using this method, and there waslittle or no dependence of accuracy on the design image. In addition,the improvement in accuracy began flattening at the 7×7 filter size anddid not improve beyond the 15×15 filter size. General purpose CPUs haveadvanced in processing power so that it is now possible to implementonline image processing using 7×7 filters without any additionalhardware. The fact that large filters are not required is advantageousfor real-time operation.

This image segmentation mechanism provides the basis of fully automated,real-time cytometry from the in-focus image stream. The advantages ofthe precision and accuracy of this method result in improved measurementfidelity, improving system throughput and resolution as explored below.

Segmentation of objects from background using least-squares-designedcontrast-enhancing filters can be used to find any cellular compartmentor set of compartments from a set of fluorescence images. Tessellationcan then be used to assign each cellular compartment to a cell. Imagesegmentation of the cell nuclei creates a single object for each cell.The nuclear masks are then fed to the tessellation algorithm to map outregions belonging to each cell and the remaining compartments areassigned to a cell based on those regions. Thus, tessellation providesan objective means to assign cellular compartments to cells.

Image segmentation of objects from background also does not separateoverlapping objects. Even in carefully controlled cell cultures, cellsmay overlap. In many cytometry applications, measurements can beimproved by cutting the overlapping objects apart. Tessellation can beused to separate overlapping objects. Once the images of the cell nucleihave been segmented using contrast enhancing filters, the nuclearpositions (e.g., centroids) can be input to the tessellation algorithmto cut the images of other overlapping cellular compartments apart andimprove measurement fidelity. Mathematical morphology techniques(erosion and dilation) and the watershed algorithm are also often usedto separate overlapping objects. For example, a cell membrane stain(e.g., DiI or DiD, Molecular Probes, Eugene Oreg.) or an amine wholecell stain (Alexa 488, Molecular Probes, Eugene Oreg.) can be used toidentify the cytoplasmic and nuclear compartments together (thecytoplasmic compartment can be determined by subtracting the nuclearcompartment). Erosion/dilation or watershed can then be used to separatethe overlapping cells. Other cellular compartments (e.g., endoplasmicreticulum, ER) or vesicles can then be assigned to the cells based onthe resulting image segmentation masks. However, very thin cellularregions of the cell that are very dim may not be segmented and may beabsent from the resulting masks. Cell compartments such as ER orvesicles that fall on the missing areas will not be assigned to a cell.Tessellation can then be used to complete the assignment. It also may beinconvenient or very difficult to add a stain to identify the cytoplasm.In the absence of a cellular mask, tessellation can be used to assigncellular compartments to cells (nuclear masks).

In this invention, once the magnified images of the cell nuclei havebeen segmented using contrast enhancing filters, tessellation of thesegmented image is utilized to precisely define nuclear positions (e.g.,centroids) in order to cut the images of other overlapping cellularcompartments apart and improve performance. Tessellation, according tothe unabridged, on line Merriam-Webster Dictionary is “a covering of aninfinite geometric plane without gaps or overlaps by congruent planefigures of one type or a few types”.

Essentially, tessellation of a magnified image in this descriptionconcerns the formation of a mosaic or a mesh of plane figures on amagnified image of cells in which each plane figure of the mosaic ormesh contains one object or compartment within a cell. Such objects orcompartments include, without limitation, nuclei, membranes, endoplasmicreticuli, Golgi, mitochondria, and other cellular regions having proteinconcentrations that are organized in some way. Further, such magnifiedimages are processed by segmentation to distinguish identical cellcomponents, such as nuclei, from the background and all other cellcomponents. The processed, segmented image is then tessellated toseparate each of the distinguished cell components from nearby andoverlapping neighbors. For example, refer to FIG. 7 in which segmentedcell nuclei are separated by a tessellation mesh of plane figureswherein each plane figure in the contains one nucleus. In this figure,for example, segmented nucleus 700 contained within the polygon 710 isseparated and therefore distinguishable from its near neighbor nucleus720 which is contained within polygon 730. Knowing the location, shape,and dimensions of each plane figure in the mosaic, the nuclei in themagnified image can be quickly accessed by traversing the mosaic in aregular fashion. The mesh in FIG. 7 is visible only in order to supportan understanding of tessellation. In practice, the mesh resulting fromtessellation can be visualized using conventional image processing,manual entry and display means. In processing of magnified images byprogrammed means, the mesh would, of course, also exist as one or morestored data structures. These data structures enable automation oftessellation by programming of the host computer 110.

The platform of FIG. 1 also includes a stable illumination source forepifluorescence. Light source stabilization can be critical in imagecytometry systems. While many other instruments (e.g.,spectrophotometers and flow cytometers) are able to overcome thisproblem by monitoring and correcting source fluctuations, imagecytometry is different because it is difficult to monitor and correctfor the 250,000 or more detectors in a CCD camera. Accordingly, thesource is constituted of a stable 100 W Hg vapor light source employingan optical fiber to scramble out the spatial variations (from arcwander), and feedback control to eliminate the remaining temporalinstabilities [S Heynen, D A Gough, and J H Price. Optically stabilizedmercury vapor short arc lamp as uv-light source for microscopy. InProceedings, International Society for Optical Engineering (SPIE),Optical Diagnostics of Biological Fluids and Advanced Techniques inAnalytical Cytology, volume 2982, pages 430-434, 1997]. The illuminationstability provided is shown in FIG. 8. The gray trace is of anunstabilized lamp and the dark trace shows the result of feedbackstabilization with our system. Both traces are envelopes showing theminimum and maximum intensity values recorded at 10 s intervals. Theintensities were measured at 30 Hz using an Imaging Technology, Inc.Series 151 image processor. This long-term stability is exceptional.With feedback, the intensity coefficient of variation was reduced byover 30 times. This intensity control will greatly improve measurementsthat depend directly on the fluorescence intensity, such as in cellfunctional assays.

The invention further utilizes image area-of-interest (AOI) datastructures and caching algorithms. One aspect of scanning cytometry isthat the inherent view of the data is a contiguous image of the slidesample, composed of tens of thousands of microscope fields and gigabytesof raw image data (typically larger than the conventional 32-bitaddressable space). Efficient access to AOIs in these data is requiredby both front-end users through the graphical interface, and byapplication programmers and assay developers. This access can only beprovided by novel image data structures and caching algorithms thatorganize disk-resident AOIs in the scanned sample, and shuffle them inand out of RAM as required. An “image table” data structure has beendeveloped [E A Hunter, W S Callaway, and J H Price. A software frameworkfor scanning cytometry. In Proceedings, International Society forOptical Engineering (SPIE), Optical Techniques in Analytical CytologyIV, San Jose, Calif., volume 3924, pages 22-28, January 2000.] toorganize and manage image data for efficient access by the mostprevalent use scenarios, including continuous (user scrollable),discontinuous (arbitrary rectangle), and database driven (queryresponse) sample browsing and queries. FIGS. 9 and 10 illustrate designsfor AOI metadata structures and associated caching algorithms. The imagetable is responsible for image area-of-interest management in a formatdesigned to facilitate efficient AOI recall based on viewport overlaygeometry. Efficient caching algorithms move disk-resident AOI data intoRAM as needed. This allows for user browsing or sequential access oftens of gigabytes of image data or more on standard PC workstations.

Automated Fluorescence Microscopy for Cell Functional Analysis in aCytoplasm-to-Nucleus NFκB Translocation Study

This section describes the results of a study to quantify thesubcellular distribution of the regulatory protein nuclear factor κB(NFκB) in response to cellular stimulation. Upon stimulation, forexample by proinflammatory cytokines, the NFκB's inhibitory subunit isphosphorylated and subsequently destroyed by proteasomes. Loss of theinhibitory subunit frees the modified protein's p65 regulatory subunitto translocate from the cytoplasm into the nucleus where it can activatedefense genes in response to the external stimulation.

Accurate and precise subcellular quantification of immunofluorescentlylabeled NFκB from microscope images provides a direct means of assessingthe ability of a compound to inhibit cellular function in the presenceof stimulation. This study examines cytoplasm-to-nucleus NFκBtranslocation in HUVEC cells in response to tumor necrosis factor α(TNFα), as an archetype of an image-based, cell functional assay.Because any other cellular compartments can be isolated and preciselyquantified by specific immunofluorescent labeling working in tandem withthe image acquisition and analysis technology described above, the speedand precision advantages demonstrated here are generally available tocell functional assays. The experiments also allow direct comparison ofthe results achieved by practice of our invention using the platform ofFIG. 1 against previously-published results using different technology[G Ding, P Fischer, R Boltz, J Schmidt, J Colaianne, A Gough, R Rubin,and D Uiller. Characterization and quantitation of NFκB nucleartranslocation induced by interleukin-1 and tumor necrosis factor-α.Journal of Biological Chemistry, 273(44):28897-28905, 1998] as measuredby fidelity and system throughput. We show comparisons of cellcompartment measurements in a substantially reduced well-populationstandard deviation, which is a function of both measurement technologyand inherent biological heterogeneity. We show how response variabilitytranslates into reductions in the number of cells required to reliablydetect fractional responses, and also how it affects inhibitor responseestimates (via a Monte Carlo simulation of a model response IC₅₀distribution). The ability to segregate cell images based onmorphological information leads to further improvements by focusingmeasurement on a homogeneously responding subpopulation. Scanning andanalysis are then validated by analyzing the response of three BIPItranslocation inhibitor compounds. Finally, system throughput isexplained and estimated scan rates are given for typical experimentalparameters.

Experimental Procedures

The cells used in these studies were primary human umbilical veinendothelial cells (HUVEC) obtained from Clonetics Corporation. The cellswere maintained in culture under standard conditions and passaged usingClonetics' proprietary EGM medium and reagent system. Since they are nottransformed, HUVECs, generally become senescent after nine or tenpassages. Cells used in these studies were from earlier passages and didnot exhibit senescence associated morphological changes.

Prior to assay, cells were transferred to 96-well plates (Packard blackViewPlate) and incubated overnight to yield cell monolayers that wereapproximately 25% confluent. Plates used to determine a statisticallyoptimal cell density for the assay were made by seeding wells at 5000,2500, and 1000 cells per well and incubating overnight. Three selectedcompounds were tested for inhibition of TNFα stimulated NFκBtranslocation in HUVEC. The three compounds and controls were laid outin the plates as shown in Tables 2 and 3. Test compounds were directlydiluted from DMSO stock solutions into medium to yield 60 mM compoundsolutions containing less than 0.7% DMSO. These solutions were seriallydiluted in medium to generate compound concentrations as low as 0.1 mM.After the medium was removed from the test plate, 100 μl aliquots ofeach compound dilution were dosed into triplicate wells. Unstimulatedand TNFα stimulated control wells received 120 and 100 μl of medium,respectively. The cells were pre-incubated for 30 minutes at 37° C.before they were stimulated by adding 20 ml of TNFα (1200 U/ml medium)to each well. The final stimulating concentration of TNFα used in thisassay was 200 U/ml. After 15 minutes incubation, the cells were fixedwith 3.7% formaldehyde and then processed by using the NFκB kit reagentsand protocol obtained from Cellomics (Pittsburgh, Pa.) to stain forcellular NFκB. In brief, cells are permeabilized, washed, incubated withrabbit anti-NFκB primary antibody for 1 hour, washed, incubated with asecondary anti-rabbit IgG antibody-Alexa Fluor 488 conjugate for 1 hour,and washed again. Nuclei were stained with either Hoechst dye includedin the secondary antibody solution or with DAPI. When used, DAPI wasadded in a final wash buffer solution at 100 to 400 ng/ml and kept inplace during storage and examination. Translocation of NFκB from thecytoplasm to the nucleus was assessed by visual examination of samplewells with a confocal microscope system.

TABLE 2 Template for 96 well NFκB control plate. Row Cells per Well ABlank B 5000 C 5000 D 2500 E 2500 F 1000 G 1000 H Blank Rows B, D, Funstimulated; rows C, E, G stimulated with 200 U/ml TNFα for 15 m.

TABLE 3 Template for 96 well BIPI inhibitor plate. Rows Rows Rows RowsRow 1-3 4-6 7-9 10-12 A 200 U/ml TNFα 200 U/ml TNFα Unstimu- Unstimu-lated lated B 0.1 μM A 0.1 μM B 0.1 μM C Blank C 0.3 μM A 0.3 μM B 0.3μM C Blank D 1 μM A 1 μM B 1 μM C Blank E 3 μM A 3 μM B 3 μM C Blank F10 μM A 10 μM B 10 μM C Blank G 30 μM A 30 μM B 30 μM C Blank H 60 μM A60 μM B 60 μM C Blank Row A controls: stimulated with 200 U/ml TNFα for15 m. (columns 1-6), unstimulated (columns 7-12). Rows B-H: triplicatewells of compounds A, B and C, at 0.1, 0.3, 1, 3, 10, 30, 60 μMconcentrations.

The high-throughput microscopy platform illustrated above in FIG. 1 wasused to automate fluorescent image acquisition and analysis of the NFκBassay. The primary system components consist of a standard fluorescentinverted microscope (Nikon TE-300); a Pentium III workstation; a controlsystem, which houses the many electronic components; and software whichcontrols the system for acquisition and analyzes the acquired images.The objective used in the scan was a Nikon CFI60 20×0.5 NA Plan Fluor.The spatially and temporally stabilized fluorescent light source is auses a 100 watt Osram HBO 103 W/2 mercury arc lamp as its light sourceand connects to the microscope at the Nikon epi-fluorescent attachment.A standard Chroma Dapi/Hoechst filter cube with a long pass emissionfilter was used for the nuclear channel. A standard FITC filter cubewith a long pass emission filter was used for the Alexa488-labeled NFκBchannel. A Cohu 640×480 pixel 10-bit progressive camera (model6612-3000; 9.9 μm² pixels) was used for image acquisition.

Measurement of Distributions or Fractional Localized Intensities (FLI)of Subcellular Compartments

Cellular substances are dynamically distributed during cellularresponses. Although there may be hundreds of thousands to millions ofdifferent cellular substances, a cellular response can be measured byspecifically labeling the subset of substances participating in theresponse. At any point in time, the combination of compartmentidentification and specifically labeled substances can be used to take asnapshot of the distributions. Image segmentation creates the imagemasks for each cellular compartment. For example, aleast-squares-designed contrast-enhancing filter is used on each of thecellular compartment images to create segmentation masks. The nuclearsegmentation is then used as a guide to separate overlappingcompartments using tessellation. These image segmentation techniqueswork best on many types of cellular images. But other segmentationtechniques can also be used to generate the segmented masks for eachcellular compartment. The measurement logic then loops through each setof pixels identified by the compartment masks and sums pixel intensitiesI(x; y). As an example, assume membrane, cytoplasmic and nuclear masks,m, c and n, respectively, with sizes N_(c), N_(n) and N_(m). Thedistributions over the compartments are defined as the fractionallocalized intensities of each compartment. The fractional localizedintensity of the cytoplasm F_(c) is

$\begin{matrix}{F_{c} = \frac{\sum\limits_{{({x,y})} \in c}{I\left( {x,y} \right)}}{\sum\limits_{{{({x,y})} \in c},m,n}{I\left( {x,y} \right)}}} & (1)\end{matrix}$

The equations are analogous for the fractional localized intensities ofthe nucleus F_(n) and membrane F_(m), and F_(c)+F_(n)+F_(m)=1. Thephysics of fluorescence image formation leads to the use of integratedintensity to quantify cellular distributions. The emission intensity atpixel location (x, y) in an image plane is

$\begin{matrix}{{I\left( {x,y} \right)} = {{\int_{0}^{z}{I_{0}Q\; ɛ\; u{z}}} = {I_{0}Q\; ɛ\; u^{\prime}z}}} & (2)\end{matrix}$

with incident (excitation) intensity I₀, quantum yield Q, extinctioncoefficient ε, local and column average fluorophore concentrations u andu′, and column thickness z. When the depth of field is greater than thecell, image formation effectively integrates the sample in z, thedirection of the optical axis. When the depth of field is smaller thanthe cell as with high NA, confocal and multiphoton optics, explicitintegration in z may more accurately represent the intensity. Assumingthis integration has already taken place either optically with largedepths of field or computationally with small depths of field, intensitymeasurements integrate in the orthogonal dimensions

$\begin{matrix}{{\sum\limits_{({x,y})}{I\left( {x,y} \right)}} = {I_{0}Q\; ɛ{\sum\limits_{({x,y})}{{u^{\prime}\left( {x,y} \right)}{z\left( {x,y} \right)}}}}} & (3)\end{matrix}$

which has units proportional to moles fluorophore. F_(c) becomes

$\begin{matrix}{F_{c} = \frac{\sum\limits_{{({x,y})} \in c}{{u^{\prime}\left( {x,y} \right)}{z\left( {x,y} \right)}}}{\begin{matrix}{{\sum\limits_{{({x,y})} \in c}{{u^{\prime}\left( {x,y} \right)}{z\left( {x,y} \right)}}} + {\frac{Q_{n}}{Q_{c}}{\sum\limits_{{({x,y})} \in n}{{u^{\prime}\left( {x,y} \right)}{z\left( {x,y} \right)}}}} +} \\{\frac{Q_{m}}{Q_{c}}{\sum\limits_{{({x,y})} \in m}{{u^{\prime}\left( {x,y} \right)}{z\left( {x,y} \right)}}}}\end{matrix}}} & (4)\end{matrix}$

with units of moles fluorophore. As before, the equations are analogousfor the fractional localized intensities of the nucleus F_(n) andmembrane F_(m). Quantum yields are potentially compartment specific(i.e., vary with pH, ion concentration and other local physicalparameters) and can be established experimentally as part of protocoldevelopment. Note that direct integration over the compartment imagesegments is preferred over average compartment intensity ratios ordifferences used in U.S. Pat. No. 5,989,835 because a ratio of areascauses a bias factor that confounds direct interpretation unlessN_(c)=N_(n)=N_(m), which can only be artificially achieved by discardinga majority of the cytoplasm and nuclear signal because typicallyN_(c),N_(n)>Nm. The cellular (or subcellular compartment) area is afunction of height even when volume is fixed. The same amount of acellular substance can be distributed over different volumes or areas.Averaging the intensities over the areas thus introduces a dependence onsomething that has nothing to do with the amount of the labeledsubstance. Note also that in equations (1) and (4) F_(c) is the fractionof the total integrated fluorescence over all of the compartments,rather than the fraction of one compartment to the other.

Generalizing further, a compartment k in an arbitrary number κ ofcompartments ζ has fractional localized intensity

$\begin{matrix}{F_{\zeta_{k}} = {{\frac{\sum\limits_{{({x,y})} \in \zeta_{k}}{I\left( {x,y} \right)}}{\sum\limits_{i = 1}^{\eta}{\sum\limits_{{({x,y})} \in \zeta_{I}}{I\left( {x,y} \right)}}}\mspace{14mu} {and}\mspace{14mu} {\sum\limits_{i = 1}^{\eta}F_{\zeta_{I}}}} = 1}} & \left( {5a} \right)\end{matrix}$

Similarly, any subset of compartments can be combined by addition toproduce multi-compartment fractional localized intensities. For example,the fractional localized intensity of compartments 1 and 3 is

F _(ζ) _(1,3) =F _(ζ) ₁ +F _(ζ) ₃   (5b)

and so on.

Results and Error in the FLIC

Fractional Localized Intensity Measurements on TNFα-induced NFκBCytoplasm-Nucleus Translocation: Immunofluorescent staining of the NFκBp65 regulatory subunit allows for easy characterization of theintercellular NFκB distribution by visual inspection of images acquiredduring scanning. The majority of fluorescence lies outside the nucleararea in unstimulated cells, while a substantial amount co-locates withthe nuclear mask subsequent to maximal stimulation with TNFα (FIG. 11).The panels of FIG. 11 show translocation of cytoplasmic NFκB p65regulatory subunit (unstimulated; left images) to cell nucleus(stimulated; right images) subsequent to stimulation of HUVEC cells withTNFα. Cell densities were 5000 cells/well (top; rows B,C), 2500cells/well (center; rows D,E) and 1000 cells/well (bottom; rows F,G). Aclear majority of Alexa488-labeled p65 resides ubiquitously in thecytoplasm before stimulation, while bound to the inhibitory subunit.Stimulation with TNFα causes phosphorylation and separation of p65 fromthe inhibitory subunit, allowing for translocation to occur. Only ≈20%of gross p65 translocates upon maximal stimulation (qualitative visualinspection of images tends to overemphasize the higher fluorophoreconcentration in the nuclei of stimulated cells, and suggest a greatertranslocation). Note how background fluorescence varies between wells.As stipulated in previous reports, the fraction of p65 to translocate onstimulus represents 20% of the gross amount contained in unstimulatedcell cytoplasm. [Ding et al.]. Effects due to background fluorescenceare corrected by estimating and subtracting the mean background imageintensity, which was found to be approximately constant fromimage-to-image, but can vary from well-to-well. Therefore, correctionswere performed on a per well basis.

Co-location of p65 immunofluorescence with the nuclear mask can occur inunstimulated cells as image formation integrates the p65-boundfluorophore emission through the three dimensional sample parallel tothe optical axis, and therefore is susceptible to contribution above orbelow the nuclear volume. Similarly, stimulated cell images carryresidual cytoplasmic p65 immunofluorescence contributions co-locatedwith the nuclear mask. Although this effect may be small, it introducesa bias in fractional localized intensity measurements. For applicationswhere this effect is not negligible, a correction procedure could bedeveloped by correlating nuclear and cytoplasm absolute integratedintensities in unstimulated cells to establish the average contributionof cytoplasm resident fluorophore to nuclear measurements as a functionof cytoplasmic intensity.

FIG. 12 illustrates image acquisition processing stages in the system ofFIG. 1. To obtain the images of FIG. 12, the nuclear channel (Hoechst33342 or DAPI) was scanned over the 10×10 field area in each well topopulate the image table with nuclear intensity images (upper left). Theupper right panel shows the associated nuclear mask images produced byfiltering and intelligent thresholding of the in-focus acquired images.Subjectively, masks are consistent with the quantitative study conductedby Price, Hunter, and Gough, which demonstrated accurate and precisenuclear masks despite significant local image variations and a widerange of nuclear intensities from cell-to-cell. Such variations canconfound simple techniques, because the accuracy and precision of theresulting masks would substantially bias localized intensitymeasurements. Using human-segmented training and testing data, thisfiltering technique has been shown to achieve accuracy >93% correctclassification of pixels and precision of less than about 1% coefficientof variation [Price, Hunter, and Gough]. This is acceptable, and perhapsa limiting performance level, since human-human agreement in manuallysegmented examples is likely to be comparable. The scan area tessellatedby nuclear centroids is shown in FIG. 12 (lower left). This tessellationstructure is then used to dynamically assemble and associate thequantitative p65 channel (Alexa488) with the individual nuclear images(lower right). Reliable alignment between channels provides the basisfor accurate fractional localized intensity measurements.

In FIG. 12, individually addressable images of Hoechst 33342- orDAPI-labeled cell nuclei were identified and extracted from the imagestream produced during scanning. These images were saved to disk usingan image table data structure (FIG. 9) that maintains the relativepositions of each nuclei (upper left). Binary nuclear mask images (upperright) show the accuracy of the nuclear segmentation processingachieved, which is critical for compartment localized measurements.Voronoi tessellation of nuclear centroids (lower left) partitioned thedata structure and scan area into polygonal regions used to correlatesubsequent channels with individual nuclei (lower right), in order tomaximize the signal from every cell.

Maximal NFκB translocation following 200 U/ml TNFα stimulation of HUVECcells was quantified by measuring fractional localized intensity in thenucleus (FLIN) (fractional localized intensity in the cytoplasm (FLIC)is equivalent in the two-compartment translocation model used in thisstudy: FLIC+FLIN=1) for every cell completely within the scan area(Table 4). At the different cell densities, FLIN sample mean, standarddeviation (σ), standard error (SE) and coefficient of variation (CV)were calculated in 12 replicate wells (Table 4). Well-to-well samplestatistics for all six density x stimulation treatments of thewithin-well sample mean FLINs are reported in Table 5, and illustratedby the constant horizontal lines in FIG. 13. FIG. 13 illustrates thepopulation statistics for TNFα-induced NFκB activation in HUVEC cells.FLIN response (μ±2σ) is plotted per well for unstimulated (lower FLIN)and maximally stimulated (increased FLIN) cells (200 U/ml TNFα) for 5000cells/well (top; rows B,C), 2500 cells/well (center; rows D,E) and 1000cells/well (bottom; rows F,G). The horizontal lines give across rowmean-of-means±2×SE.

TABLE 4 Well-wise and row-pooled FLIN statistics. CVs defined as (stddev)/mean Well n mean σ SE CV  1B 964 0.147846 0.048122 0.0015500.325484  2B 1132 0.130526 0.042521 0.001264 0.325768  3B 1309 0.1344700.043876 0.001213 0.326291  4B 1263 0.132880 0.043633 0.001228 0.328090 5B 1089 0.137363 0.045115 0.031367 0.328437  6B 1297 0.138022 0.0451320.001253 0.326990  7B 1251 0.138558 0.050421 0.001426 0.363893  8B 11630.134308 0.046321 0.001364 0.344884  9B 1283 0.132846 0.046153 0.0012890.347419 10B 1127 0.132862 0.044833 0.001335 0.337442 11B 994 0.1261780.040244 0.001276 0.318942 12B 1234 0.129016 0.043913 0.001250 0.340366Row 14096 0.1345 0.0454 0.0004 0.3374   B  1C 1123 0.290729 0.0757130.002259 0.260423  2C 1226 0.286024 0.076772 0.002193 0.268410  3C 14430.318096 0.085135 0.002241 0.267640  4C 1265 0.346417 0.077291 0.0021730.223114  5C 1410 0.342634 0.081898 0.002181 0.239025  6C 1149 0.3388050.081211 0.002396 0.239627  7C 1367 0.308078 0.082221 0.002224 0.266019 8C 1330 0.305640 0.081423 0.002233 0.266401  9C 1260 0.314321 0.0804580.002267 0.255975 10C 1311 0.322215 0.080099 0.002212 0.248569 11C 11080.328182 0.083463 0.002507 0.254319 12C 1044 0.336819 0.082491 0.0025530.244912 Row 15036 0.3199 0.0829 0.0007 0.2592   C  1D 192 0.1653600.055521 0.004007 0.335721  2D 228 0.169188 0.066621 0.004412 0.393770 3D 273 0.177909 0.062334 0.003773 0.350374  4D 330 0.167759 0.0637490.003509 0.380006  5D 356 0.219419 0.086005 0.004558 0.391968  6D 3670.211853 0.062027 0.004262 0.387188  7D 349 0.184329 0.069234 0.0037060.375601  8D 319 0.183501 0.068811 0.003853 0.374990  9D 311 0.2193870.092681 0.005255 0.422455 10D 382 0.150226 0.052488 0.002686 0.34939011D 354 0.158436 0.058115 0.003088 0.366802 12D 277 0.149743 0.0513060.003083 0.342629 Row 3738 0.1809 0.0735 0.0006 0.4062   D  1E 4880.350998 0.088568 0.004009 0.252331  2E 533 0.352654 0.076766 0.0033250.217679  3E 591 0.350479 0.074648 0.003071 0.212983  4E 562 0.3710060.077932 0.003287 0.210056  5E 615 0.368168 0.077079 0.003108 0.209359 6E 589 0.361919 0.080731 0.003326 0.223064  7E 594 0.361462 0.0773520.003174 0.213998  8E 542 0.369526 0.081884 0.003517 0.221593  9E 6220.359408 0.089347 0.003473 0.248594 10E 644 0.352375 0.079729 0.0031420.226263 11E 534 0.367225 0.082335 0.003563 0.224207 12E 396 0.3740300.091069 0.004576 0.243482 Row 6750 0.3613 0.0816 0.0010 0.2259   E  1F199 0.171691 0.052991 0.003756 0.308639  2F 109 0.154424 0.0750780.007191 0.486179  3F 204 0.147685 0.056180 0.003933 0.380404  4F 1820.149932 0.049610 0.003677 0.330882  5F 177 0.153247 0.058143 0.0043700.379411  6F 174 0.155897 0.051656 0.003916 0.331348  7F 200 0.1516080.056382 0.003987 0.371893  8F 172 0.183963 0.068814 0.005247 0.374063 9F 163 0.146976 0.056036 0.004389 0.381262 10F 171 0.146951 0.0499640.003821 0.340006 11F 176 0.155161 0.051188 0.003858 0.329902 12F 1600.158965 0.049600 0.003921 0.311977 Row 2087 0.1564 0.0570 0.0012 0.3642  F  1G 279 0.336876 0.086114 0.005155 0.255624  2G 332 0.4064670.077388 0.004247 0.190391  3G 306 0.340058 0.078636 0.004495 0.231235 4G 296 0.348505 0.068383 0.003975 0.195218  5G 310 0.359116 0.0703160.003994 0.195803  6G 339 0.346658 0.077604 0.004215 0.223862  7G 3290.354841 0.075766 0.004177 0.213521  8G 290 0.349761 0.074941 0.0044010.214263  9G 281 0.362283 0.075045 0.004477 0.207145 10G 292 0.3445760.077994 0.004564 0.226347 11G 317 0.341355 0.076450 0.004294 0.22396112G 235 0.347623 0.069314 0.004522 0.199394 Row 3606 0.3537 0.07790.0013 0.2203   G

We measured the occurrence of 18.18%-19.68% (18.80% average)translocation of labeled NFκB, which we calculated by averaging the 12row replicate well FLIN sample means per row, and differencingstimulated and unstimulated averages at each cellular density.Heterogeneity of cellular response to TNFα stimulation is apparent byvisual inspection of acquired images and summarized in the ±2σconfidence interval widths in FIG. 13 (the larger interval around eachwell mean). Stimulated and unstimulated well populations havesignificant overlap in response. However, the average response of 18.80%is large compared to the FLIN standard errors which are calculatedwithin wells in Table 4 and shown graphically in FIG. 13 (smallerinterval about each well mean). For example, the average SE in row E,3.46×10⁻³, gives a 4 SE empirical response resolution of 0.0139 FLINincrements, or about 7.4% increments over the 18.80% signal dynamicrange.

By analyzing the spread of well-average FLIN measurements fromwell-to-well, we can assess well-to-well variability and repeatabilityof the experiment. The row aggregate FLIN sample means, calculated byaveraging the 12 replicate well sample means in each row, have astandard error of 1.6×10⁻³-7.20×10⁻³ (CV of 2.34%-13.93%) (Table 5).This well-to-well sample mean variability is visualized by thehorizontal aggregate mean and ±2SE lines in FIG. 13. Deviations of theindividual well responses from this region show the effect of protocolperturbations in individual well environments. Although statisticallysignificant differences exist, the physical or scientific significancein terms of effect on measured translocation is typically small in thisexperiment, indicating the possibility of limiting the number ofreplicate wells required in an optimized and controlled protocol. As aquantitative example, a trimmed between—well CV for row D in Table 4reduces variability from 13.93% to 4.76%. Note that row D had thehighest CV and that the numbers were trimmed to the six median wells(using the assumption that those wells with the most human preparationerror were at the tails of the distribution). Wellplate replicate layoutto mitigate the effect of well-to-well variability must be determinedduring specific assay development.

TABLE 5 Well-to-well statistics of well-average FLIN values (n = 12).Row mean σ SE CV B 0.1346 0.0055 0.0016 0.0412 C 0.3199 0.0198 0.00570.0618 D 0.1798 0.0250 0.0072 0.1393 E 0.3616 0.0085 0.0024 0.0234 F0.1564 0.0110 0.0032 0.0704 G 0.3532 0.0184 0.0053 0.0521 Mean, standarddeviation, standard error and coefficient of variation of the FLIN wellsample means across each row.

Automatic Determination of Cells/Well Required to Reach a MinimumSignificant Response

Theory of Dose Response Resolution: To underpin the empirical doseresponse resolution theoretically and to properly stage a baselinecomparison of the data with results produced using a differenttechnology, an objective definition of the meaning of responseresolution and its dependence on measurement fidelity, sample size anderror controls is needed. Inhibitory responses are estimated bynonlinear regression from a collection of experimental well-averagemeasurements distributed across a range of inhibitor concentrations.Sources of variability in curve estimates include high populationvariability (heterogeneity) in response, well-to-well variability andinadequate sampling of inhibitor concentration. These factors arenaturally dependent; measurement precision places limitations onmeasurable differences in response. Experiment design involvesoptimizing jointly over well datapoint replicates and individual wellmeasurement quality (a function of the number of cells per well).

A direct measure of system resolution is the minimum statisticallysignificant response that can be reliably measured. An equivalentmeasure of system performance is, for a specific minimum significantresponse implicated by choice of inhibitory response model, the numberof cell measurements required at each inhibitor concentration toreliably estimate the response model parameters. For two populations ofresponding cells, we want to determine the minimum difference in meanFLIN that can be reliably detected as a measure of system resolution.For this, we define a two-sample hypothesis test

H ₀:μ₁−μ₂=0  (6)

H _(α):μ₁−μ₂>0  (7)

where μ₁ and μ₂ are the true (unobserved) mean FLIN responses, and H₀and H_(α) are the null and alternative hypotheses. This is anupper-tailed test with decision rule z>z_(α), where z=(x₁−x₂)/(σ√2/n) isthe unit variance normalized test statistic for sample means x₁ and x₂,and z_(α) is the threshold for rejection of an α-level test,

$\begin{matrix}{\alpha = {{\Pr \left( {z > z_{\alpha}} \middle| H_{0} \right)} = {{1 - {\int_{- \infty}^{z_{\alpha}}{\frac{1}{\sqrt{2\pi}}^{- \frac{x^{2}}{2}}{x}}}} = {1 - {\Phi \left( z_{\alpha} \right)}}}}} & (8)\end{matrix}$

The type I error probability (α) is the likelihood that two identicallyinhibited samples will, by natural variations, show a minimallysignificant response difference in measurements). The probability ofType I errors is controlled by fixing its value in advance, e.g. α=0.05,as determined by assay requirements. Similarly, the type II errorprobability (β) is the likelihood that two samples with minimallydifferent responses will not show a significant difference inmeasurements, and it is similarly determined and fixed by assayrequirements. To relate these assay parameters and derive a measure ofprecision, we express the probability β mathematically, assuming anabsolute minimum significant response (MSR) Δμ

$\begin{matrix}{\mspace{79mu} {{\Pr \left( {\left. {z < z_{\alpha}} \middle| {\mu_{1} - \mu_{2}} \right. = {\Delta \; \mu}} \right)} = {\beta ({\Delta\mu})}}} & (9) \\{{\Pr\left( {\left. {\frac{x_{1} - x_{2} - {\Delta\mu}}{\sigma \sqrt{2/n}} < {z_{\alpha} - \frac{\Delta\mu}{\sigma \sqrt{2/n}}}} \middle| {\mu_{1} - \mu_{2}} \right. = {\Delta\mu}} \right)} = {\beta \left( {{\Delta\mu},\sigma,n} \right)}} & (10) \\{\mspace{79mu} {{\Phi\left( {z_{\alpha} - \frac{\Delta\mu}{\sigma \sqrt{2/n}}} \right)} = {\beta \left( {{\Delta\mu},\sigma,n} \right)}}} & (11)\end{matrix}$

This expresses β as a function of Δμ, σ (FLIN population standarddeviation) and n (sample size). The assumption that both populationsshare the same σ is reasonable because in the limit, the minimumdetectable difference approaches the same population (formulae for thecase σ₁≠σ₂ are easy to derive, however, when we need an estimate ofsample size requirements for detecting very distinct responses). Byspecifying β as we did for α, e.g. α=β=0.05, we control type II errorprobability and fix z_(β)

$\begin{matrix}{{\Phi \left( z_{\beta} \right)} = \beta} & (12) \\{z_{\beta} = {z_{\alpha} - \frac{\Delta\mu}{\sigma \sqrt{2/n}}}} & (13) \\{{\Delta\mu} = {\left( {z_{\alpha} - z_{\beta}} \right)\sigma \sqrt{2/n}}} & (14)\end{matrix}$

MSR is expressed as a fraction of the dynamic range of the assay orexperiment

$\begin{matrix}{{MSR} = {\frac{\Delta\mu}{\mu_{m\; {ax}} - \mu_{m\; i\; n}} = {\frac{\left( {z_{\alpha} - z_{\beta}} \right)\sigma \sqrt{2/n}}{\mu_{m\; {ax}} - \mu_{m\; i\; n}}.}}} & (15)\end{matrix}$

For specified MSR and assay parameters, the minimum corresponding samplesize is

$\begin{matrix}{n = {2\left( \frac{\left( {z_{\alpha} - z_{\beta}} \right)\sigma}{{\Delta\mu}\; {MSR}} \right)^{2}}} & (16)\end{matrix}$

To gain an understanding of this precision measure, let α=β=0.05, thenMSR is the minimum significant dose response relative to the assaydynamic range such that

-   -   1. There is a 95% chance that when there is no true difference        between population mean responses, z≦z_(α) or        x₁−x₂≦z_(α)σ√2/n=Δμ/2 (no significant difference measured        between sample means), and    -   2. There is a 95% chance that when μ₁−μ₂=Δμ>0, x₁−x₂>Δμ/2        (significant difference measured between sample means)

Therefore, specification of MSR and protocol parameters allows one tomake objective guarantees about the minimum required sample size tocontrol the variability of dose response point measurements so that theyare unlikely to overlap (according to our specification of α, β). Afamily of n(MSR) curves are plotted in FIG. 14. The average population σin row E is σ=8.1453×10⁻², which would allow a 7.4% MSR for n≈742.

The control of type I and II error probabilities in this scheme hasdifferent implications than most traditional hypothesis testing byvirtue of its use as a measure of experimental precision. A type I errorin an inhibitor assay indicates that the replicate measurements producedresponses outside the experimentally set control interval. A type IIerror indicates that nonreplicate measurements (different inhibitorconcentrations) produced measured responses within the same controlinterval. Nonlinear regression to an inhibitory response model will notdiscriminate the different meanings of these errors; only that they bothproduce datapoint scatter and complicate parameter estimation.Therefore, for experiments involving inhibitory response, it issuggested to constrain α=β.

FIG. 14 is a curve family showing theoretical dose response resolutionfor σ={0.02,0.04,0.06,0.08,0.1,0.12,0.14} (ascending, α=0.05 (Type Ierror probability), β=0.05 (Type II error probability)). In this figure,the number of cells required (y axis) is plotted against the fraction oftotal NFκB translocation that is resolved (x axis). Measurements madeaccording to this invention fall between σ=0.06 and σ=0.08 (betweenthird and fourth curves from the bottom) as reported in Table 4. Theserelationships allow for easy visualization of the effect of improved ordegraded measurement fidelity and resolution on cell numberrequirements.

Sample Size Requirements: Ding et al. [G Ding, P Fischer, R Boltz, ISchmidt, I Colaianne, A Gough, R Rubin, and D Hiller. Characterizationand quantitation of NFκB nuclear translocation induced by interleukin-1and tumor necrosis factor-α. Journal of Biological Chemistry,273(44):28897-28905, 1998] report sample size requirements for similarNFκB translocation experiments carried out on an imaging platform. Todemonstrate the superior cell measurement fidelity of the technologyreported here, Table 6 compares the results published by Ding et al.against similar calculations based on the raw measurements data in Table4 and the derivation in the previous subsection. Table 6 providescomparative estimates of the number of cellular measurements needed tomeasure a given level of stimulation response (from 2%-100% of the totaltranslocation amplitude of about 19% FLIN change) as a function of NFκBfluorophore translocation at 95% and 99% type I confidences(α=0.05,0.01; β=0.20). The parameters of the experiment were duplicatedhere to provide a direct comparison to on cellular measurement fidelity.Additional response level 2%-10% and an additional column (α=β=0.001)were added to further demonstrate results according to our invention. Inthe general population, the same number of cells required to measure a100% response (19% translocation event) using the technology reported inDing et al. allowed a measurement of about 25% response using theprinciples of our invention. Reasonable dose response curve estimationmay require a 20% resolution or better, which at α=0.05, β=0.20 wouldrequire about 12× fewer measurements to achieve using the principles ofour invention.

TABLE 6 Comparison of the minimum number of cells required to measurepercent of the total (maximum stimulation) NFκB translocation (19% FLINchange) obtained with the invention with those reported in Ding et al.,at controlled error probabilities (type I, type II). % Ding et al.Invention Ding et al. Invention Invention Maximum (0.05, (0.05, (0.01,(0.01, (0.001, Response 0.20) 0.20) 0.20) 0.20) 0.001)  2% Not given3970 Not given 6445 24528  5% Not given 640 Not given 1039 3952 10% Notgiven 162 Not given 263 1000 20% >500 42 >500 68 256 30% 388 <30 >500 31117 35% 220 <30 331 <30 87 45% 141 <30 213 <30 54 55% 99 <30 149 <30 3765% 73 <30 111 <30 <30 70% 57 <30 86 <30 <30 80% 45 <30 68 <30 <30 90%37 <30 56 <30 <30 100%  31 <30 47 <30 <30 The response column gives thepercent response of the full dynamic range of the measurement (18.8%).Our calculations assume well sample mean measurements Gaussiandistributed, justified by Central Limit Theorem, for n ≧ 30. For n < 30,student's t distribution is appropriate. However, for simplicity, aminimum of 30 cells is always measured. Our calculations assume Δμ =0.1818 and σ linearly interpolated between 0.0650 (unstimulated) and0.0800 (stimulated).

Monte Carlo Response Curve Estimates: To quantify the effect ofmeasurement precision directly on estimated inhibitor responseparameters, a Monte Carlo simulation was carried out using the modelresponse shown in FIG. 15 which shows a theoretical dose response curvewith IC₅₀=8.895×10⁻⁵, and μ_(max), μ_(min) and σ specified by themeasurement procedure to be simulated. For FIG. 15, μ_(max)=1.0, andμ_(min)=0.0. The circled points indicate the 11 inhibitor concentrationsused in the simulation. These eleven sites were chosen logarithmicallyacross the inhibitor concentration range to which Gaussian random noisewith deviation σ/√n was added to create 500 simulated experimentaloutcomes. FIGS. 16 and 17 illustrate the simulated datasets, eachnonlinearly regressed to the model in FIG. 15 using the Gauss-Newtonprocedure [See, for example, W H Press, S A Teukolsky, W T Vetterling,and B P Flannery. Numerical Recipes in C, Second Edition. CambridgeUniversity Press, New York, N.Y., 1992.].

The panels of FIG. 16 comparatively illustrate worst case simulated doseresponse experiments achieved by way of empirical curve IC₅₀distribution for well-well mean FLIN response. The top three panels showestimated IC₅₀ scatterplots with true value. The middle three panels areIC₅₀ frequency histograms (compare with the 5 and 95-percentile valuesin Table 7). The bottom panels illustrate a sample of 20 estimated doseresponse curves from the total simulation of 500. The leftmost threepanels show a 90% IC₅₀ confidence interval in the worst case for ourinvention with full population, with σ=0.093. However, our worst case isapproximately 3.62× better (narrower) than the rightmost three panels(Ding et al.). Results for our G₁ subpopulation (center; worst caseσ=0.066) exploit the homogeneity of this subpopulation to give a 90%IC₅₀ confidence interval that is 5.07× better. Experimental parametersfor the simulation are given in Table 7.

In FIG. 17 best case simulated dose response experiments illustrateempirical curve IC₅₀ distribution for well-well mean FLIN response.Shown are estimated IC₅₀ scatterplots with true value (the top threepanels), IC₅₀ frequency histograms (the middle three panels; comparewith 5 and 95-percentile values in Table 7), and a sample of 20estimated dose response curves from the total simulation of 500 (in thebottom three panels). The 90% IC₅₀ confidence interval for the left handside (simulation of results produced by the invention; best case fullpopulation, σ=0.040) is approximately 8.81× better (narrower) than theright (simulation of results produced according to Ding et al.). Resultsfor the invention's G₁ subpopulation (center; best case, σ=0.021)exploit the homogeneity of this subpopulation to give a 90% IC₅₀confidence interval that is 16.3× better. Experimental parameters forthe simulation are given in Table 7.

Standard errors for response point estimates were 0.0093 (worst caseTable 4 standard deviation), 0.004 (best case Table 4 standarddeviation) and 8.354 (calculated from the cell numbers for statisticalsignificance reported in Ding et al. and replicated in Table 6); bothstandard errors assume n=100. Table 7 reports the 90% Monte Carloconfidence interval widths for all three outcomes, 5.8585×10⁻⁵,2.4063×10⁻⁵ and 2.1202×10⁻⁴, showing 3.62× (worst case) and 8.81× (bestcase) stronger confidence in the IC₅₀ estimates obtained for ourinvention. Further simulations shown in FIGS. 18 and 19 demonstrate thesensitivity of IC₅₀ estimate precision on the sampling of experimentalconcentrations.

For example, in FIG. 18, Monte Carlo simulations of dose responseregression reveal that IC₅₀ estimates becomes more sensitive tocell-measurement precision as the range of experimental concentrationsdiminishes. The panels in the left column show the model sigmoids usedin the experiments, ranging from 11 concentrations (top; circles) to 3concentrations (bottom; circles). 20 of 500 runs are graphically plottedfor the invention (middle column; σ=0.0665, Δμ=0.18), and Ding et al.(right column; σ=83.545, Δμ=56).

Refer further to FIGS. 19 a and 19 b, which show IC₅₀ Monte Carloconfidence interval width for Ding et al. (circles; σ=83.545, Δμ=56) andthe invention (plusses; σ=0.0665, Δμ=0.18) versus number of sampleconcentrations (3,5,7,9,11). FIG. 19 b is zoom of FIG. 19 a. Manifestly,measurement precision becomes more important as sampling and coverage ofthe relevant concentration ranges decreases. The elbow in the graphcorresponds to a regression breakdown, showing robust regression resultsmaintained by the measurement precision of the invention.

TABLE 7 Monte Carlo simulated dose response parameters and estimationresults for 500 runs at n = 100 cellular measurements, showing ΔIC₅₀ 90%Monte Carlo confidence interval variability under different measurementprecisions. std dev range CV IC₅₀ 5% IC₅₀ 95% ΔIC₅₀ 0.093 [.15, .33] 51% 6.2930 × 10⁻⁵ 1.2151 × 10⁻⁴ 5.8585 × 10⁻⁵ 0.040 [.15, .33] 22.2%7.8407 × 10⁻⁵ 1.0247 × 10⁻⁴ 2.4063 × 10⁻⁵ 0.066 [.15, .33] 35.1% 7.0982× 10⁻⁵ 1.1109 × 10⁻⁴ 4.1808 × 10⁻⁵ 0.021 [.15, .33] 11.7% 8.2296 × 10⁻⁵9.5332 × 10⁻⁵ 1.3036 × 10⁻⁵ 83.545  [0, 55]  152% 3.3819 × 10⁻⁵ 2.4584 ×10⁻⁴ 2.1202 × 10⁻⁴ Well-to-well mean FLIN response standard errors:0.0093 (invention; worst case full cellular population), 0.0040(invention; best case full cellular population), 0.0066 (invention;worst case G₁ subpopulation), 0.0021 (invention; best case G₁subpopulation), and 8.3545 (Ding et al.).

Homogeneous Subpopulation Analysis

To illustrate the ability of the invention to estimate quantities ofinterest conditionally on specific morphologic subpopulations, anempirically determined classification scheme was developed to isolate G₁cells from the S, M, G₂ cells and fluorescent artifacts. A six-ruleclassifier (Table 9) was developed from the set of standard cellularmeasurements (Table 8) taken in the Hoechst 33342 or DAPI nuclearchannel. Analysis of translocation experiment data was repeated for thisG₁ subpopulation at 2500 cells/well (rows D,E) and is reported in Tables10 (raw well data), 11 (statistics of well means) and FIG. 20 in whichthe translocation response of the G₁ subpopulation is plotted. The rowaggregate (average of all well average subpopulation FLIN responses) are0.1418 (unstimulated) and 0.3317 (stimulated).

The improved homogeneity of cell response is clearly seen in reduced ±2σpopulation intervals (compare FIG. 13, middle, with FIG. 20). Thesubpopulation defined by Table 9 responds about twice as uniformly asthe full cellular population, while the mean response is nearlyidentical. Well-to-well variability appears consistent with the fullpopulation analysis.

Table 12 reports sample size requirements for the NFκB translocationexperiment assuming Measurements are restricted to the G₁ subpopulation.The to improvement over the full population requirements at 20% responseresolution (giving the minimum of 5 samples used in a dose responseregression), is 42 cells to <30 cells (29% reduction) at α=0.05, β=0.20;68 to <30 (56% reduction) at α=0.01, β=0.20; and 256 to 47 (82%reduction) at α=0.001, β=0.001. The subpopulation analysis reduced cellnumber requirements by the 94% over data reported by Ding et al. (α=0.05and 0.01 for 20% response case).

The Monte Carlo simulations were also repeated using the G₁subpopulation statistics. Experimental parameters and numerical IC₅₀ 90%confidence interval widths given in Table 7 show a 28.6%-45.8% reductionin interval width over full population analysis for the same number ofcells n=100. Comparison to the simulations based on Ding et al. (asdescribed above) show a 5.07-fold improvement in the worst case,visualized in FIG. 16 (compare middle and right columns), and a16.3-fold improvement in the best case, visualized in FIG. 17.

These results indicate that improved response can be exploited byreducing the required scan area and keeping a specified resolution(optimizing for system throughput), or fixing the scan area and thusthroughput, but discarding cellular measurements in an intelligentmanner to select homogenous response populations and improving responseresolution. These advantages will only be available when the fraction ofcells that define the subpopulation exceeds a threshold defined by theresponse resolution equations. For example, a scan area can be reducedfrom the minimum required to measure the full population n_(F) whenever

$\begin{matrix}{{n_{s} < {n_{Ff}\mspace{14mu} {or}\mspace{11mu} \frac{n_{s}}{n_{F}}}} = {\frac{\sigma_{s}^{2}}{\sigma_{F}^{2}} < f}} & (17)\end{matrix}$

where f is the fraction of cells comprising the homogenoussubpopulation. When this is true, the scan area can be reduced by somenumber of fields dependent on cellular density and f. Thus,subpopulation analysis can directly affect system throughput whenever ahomogenous enough and large enough subpopulation exists. Similarly,minimum significant response is improved whenever

MSR_(S)−MSR_(F)>0 or σ_(S)−σ_(F) √{square root over (f)}>0  (18)

These rules can be validated by examining the curves in FIG. 14.

TABLE 8 Standard cellular measurements measurement definition  1 serialnumber

 2 x coordinate of nuclear bounding box center x_(b)  3 y coordinate ofnuclear bounding box center y_(b)  4 x size of nuclear bounding boxΔx_(b)  5 y size of nuclear bounding box Δy_(b)  6 nuclear area

 = # nuclear pixels  7 root nuclear area  8 integrated nuclear intensity$I_{n} = {\sum\limits_{n}{I\left( {x,y} \right)}}$  9 average nuclearintensity

 = I_(n)/

10 minimum nuclear intensity  p₀ = min[I_(n)(x, y)] 11 5th percentile ofnuclear intensity  p₅ = 5%[I_(n)(x, y)] 12 25th percentile of nuclearintensity  p₂₅ = 25%[I_(n)(x, y)] 13 50th percentile of nuclearintensity  p₅₀ = 50%[I_(n)(x, y)] 14 75th percentile of nuclearintensity  p₇₅ = 75%[I_(n)(x, y)] 15 95th percentile of nuclearintensity  p₉₅ = 95%[I_(n)(x, y)] 16 maximum nuclear intensity p₁₀₀ =max[I_(n)(x, y)] 17 interquartile range of nuclear intensity I

 = p₇₅ − p₂₅ 18 variance of nuclear intensity 19 standard deviation ofnuclear intensity

20 third central moment of nuclear intensity

 = (

 − 1)⁻¹Σ_(n)(I(x, y) −

21 3^(rd) root of 3^(rd) central moment of nuclear intensity m_(n,3)^(1/3) 22 fourth central moment of nuclear intensity

 = (

 − 1)⁻¹Σ_(n)(I(x, y) −

23 4^(th) root of 4^(th) central moment of nuclear intensity m_(n,4)^(1/4) 24 nuclear perimeter

 = # nuclear pixels (8-connected) 25 nuclear wiggle w =

/a_(n) 26 normalized nuclear wiggle w = I_(n)/{square root over (a_(n))}27 channel-j integrated intensity

 = Σ

(x, y) 28 channel-j integrated nuclear intensity ? = ?(x, y)?indicates text missing or illegible when filed 29 channel-j integratedcytoplasmic intensity ? = ?(x, y)?indicates text missing or illegible when filed 30 channel-j FLIN

/(

 +

)

indicates data missing or illegible when filed

TABLE 9 Six-rule definition for G₁ subpopulation based on cellularmeasurements defined in Table 8 rule effect (property of cells excluded)1 22.73 ≦ {square root over (a_(n))} ≦ 31.85 nuclear area is atypical ofresting cell nuclei 2 10⁵ ≦ I_(j) ≦ 4 × 10⁵ little or no NFκBfluorophore or large, bright fluorescent artifacts 3 m_(n, 3) < 0dividing, have unusual nuclear shapes or are dying 4 .11 ≦ w ≦ .2 two ormore aggregate nuclei or irregularly shaped nuclei 5 10 ≦ IQ_(n) ≦ 60irregular nuclear textures 6 w < 4.8 irregularly shaped nuclei

TABLE 10 Well-wise and row pooled FLIN statistics of G₁ subpopulation,as defined by the classifier in Tables 8 and 9. CVs defined as (stddev)/mean. These results show significant reduction in cell responsevariability, indicating a much more homogenously respondingsubpopulation amongst the total cell population in each well. See FIG.20 for a visualization of this data. Well n mean std dev std err CV D153 0.139494 0.026964 0.003704 0.193301 D2 79 0.135942 0.021481 0.0024170.158018 D3 85 0.142980 0.029554 0.003206 0.206699 D4 128 0.1373860.027846 0.002461 0.202687 D5 70 0.154410 0.027970 0.003343 0.181141 D675 0.153748 0.027003 0.003118 0.175630 D7 88 0.142782 0.025264 0.0026930.176942 D8 102 0.145129 0.029421 0.002913 0.202720 D9 72 0.1537870.025640 0.003022 0.166723  D10 114 0.127378 0.023832 0.002232 0.187094 D11 106 0.136946 0.026416 0.002566 0.192890  D12 107 0.131305 0.0257260.002487 0.195925 Row D 1029 0.1406 0.0277 0.0008 0.1969 E1 204 0.3281620.064145 0.004491 0.195469 E2 236 0.330019 0.058368 0.003799 0.176861 E3257 0.329866 0.056664 0.003535 0.171778 E4 230 0.350588 0.0666550.004395 0.190124 E5 269 0.349991 0.060755 0.003704 0.173590 E6 2550.324779 0.055980 0.003506 0.172363 E7 241 0.337747 0.058013 0.0037370.171766 E8 227 0.327979 0.060758 0.004033 0.185250 E9 290 0.3221160.061913 0.003636 0.192207  E10 278 0.325874 0.058762 0.003524 0.180322 E11 204 0.331738 0.064243 0.004498 0.193656  E12 162 0.320942 0.0472530.003713 0.147232 Row E 2853 0.3318 0.0606 0.0011 0.1821

TABLE 11 Well-well statistics of well subpopulation average FLIN (n =12), as defined by the classifier in Tables 8 and 9. Row mean std devstd err CV D 0.1418 0.0088 0.0026 0.0624 E 0.3317 0.0097 0.0028 0.0295Mean, standard deviation, standard error and coefficient of variation ofthe FLIN well sample means across each row. Although standard error isincreased due to the reduced number of participating cells in thesubpopulation average, well-well variability is consistent with the fullpopulation data reported in Table 5. See FIG. 20 for a visualization ofthis data.

TABLE 12 Minimum number of G₁ subpopulation cells required to resolvepercent of the total NFκB translocation at maximum stimulation, reportedhere and in [6], at controlled error probabilities (type I, type II). %Ding et al. Invention Ding et al. Invention Invention Maximum (0.05,(0.05, (0.01, (0.01, (0.001, Response 0.20) 0.20) 0.20) 0.20) 0.001)  2%Not given 624 Not given 1013 3855  5% Not given 103 Not given 168 63710% Not given <30 Not given 44 168 20% >500 <30 >500 <30 47 30% 388<30 >500 <30 <30 35% 220 <30 331 <30 <30 45% 141 <30 213 <30 <30 55% 99<30 149 <30 <30 65% 73 <30 111 <30 <30 70% 57 <30 86 <30 <30 80% 45 <3068 <30 <30 90% 37 <30 56 <30 <30 100%  31 <30 47 <30 <30 Ourcalculations assume Δμ = 0.1818 and σ linearly interpolated between0.0267 (unstimulated) and 0.0542 (stimulated). See Table 6 caption foradditional assumptions.

FIGS. 21 a and 21 b illustrate inhibition of NFκB nuclear translocationby BIPI compound A. In these figures, Fractional Localized Intensity inthe Nucleus (FLIN) is plotted versus inhibitor concentration for fullcell population found in the 10×10-field scan area (FIG. 21 a), and in a100-cell population (FIG. 21 b). Cell preparation is as described above.FIG. 21 a illustrates total cell population analysis with inhibition ofNFκB nuclear translocation in HUVEC cells where Ki(Cmpd_A)=6.95E-05,Rmax=3.86E-01, and Rmin=0.212. FIG. 21 b illustrates 100-cell populationanalysis with inhibition of NFκB nuclear translocation in HUVEC cellswhere Ki(Cmpd_A)=8.09E-05, Rmax=4.21E-01, and Rmin=0.212.

FIGS. 22 a and 22 b illustrate inhibition of NFκB nuclear translocationby BIPI compound B. In these figures, Fractional Localized Intensity inthe Nucleus (FLIN) is plotted versus inhibitor concentration for fullcell population found in the 10×10-field scan area (FIG. 22 a), and in a100-cell population (FIG. 22 b). Cell preparation is as described above.FIG. 22 a illustrates total cell population analysis with inhibition ofNFκB nuclear translocation in HUVEC cells where Ki(Cmpd_B=5.08E-05,Rmax=3.77E-01, and Rmin=0.212. FIG. 22 b illustrates 100-Cell populationanalysis with inhibition of NFκB nuclear translocation in HUVEC cellswhere Ki(Cmpd_B)=6.65E-05, Rmax=4.12E-01, and Rmin=0.212. FIGS. 23 a and23 b illustrate inhibition of NFκB nuclear translocation by BIN compoundC. Fractional Localized Intensity in the Nucleus (FLIN) is plottedversus inhibitor concentration for full cell population found in the10×10-field scan area (FIG. 23 a), and in a 100-cell population (FIG. 23b). Cell preparation is as described above. FIG. 23 a illustrates totalcell population analysis with inhibition of NFκB nuclear translocationin HUVEC cells where Ki(Cmpd_C=1.01E-06, Rmax=3.99E-01, andRmin=0.212E-01. FIG. 23 b illustrates 100-cell population analysis withinhibition of NFκB nuclear translocation in HUVEC cells whereKi(Cmpd_C)=1.36E-06, Rmax=4.20E-01, and Rmin=0.212.

FIG. 24 illustrates inhibition of NFκB nuclear translocation by BINcompound B in the G₁ cell subpopulation defined by Tables 8 and 9 withFractional Localized Intensity in the Nucleus (FLIN) plotted versusinhibitor concentration for full cell population found in the10×10-field scan area. Cell preparation is as described above. In thisfigure, G1 cell population analysis is shown with inhibition of NFκBnuclear translocation in HUVEC cells where Ki(Cmpd_B)=3.70E-05,Rmax=3.82E-01, and Rmin=0.212.

TABLE 13 Inhibitor response parameter estimates and standard errors,given by nonlinear regression analysis. Compound R_(max) (FLIN) (SE)R_(min) (FLIN) (SE) IC₅₀ (FLIN) (SE) A 0.386 (8.17 × 10⁻³) 0.212* 6.95 ×10⁻⁵ (2.2 × 10⁻⁵) A (n = 100) 0.421 (6.70 × 10⁻³) 0.212* 8.09 × 10⁻⁵(1.8 × 10⁻⁵) B 0.377 (3.63 × 10⁻³) 0.212* 5.08 × 10⁻⁵ (6.8 × 10⁻⁶) B (n= 100) 0.412 (5.85 × 10⁻³) 0.212* 6.65 × 10⁻⁵ (1.3 × 10⁻⁵) B (G₁) 0.382(3.28 × 10⁻³) 0.212* 3.70 × 10⁻⁵ (4.1 × 10⁻⁶) C 0.399 (6.90 × 10⁻³)0.212 (3.78 × 10⁻³) 1.02 × 10⁻⁶ (1.6 × 10⁻⁷) C (n = 100) 0.424 (8.96 ×10⁻³) 0.212* 1.15 × 10⁻⁶ (2.2 × 10⁻⁷) Unity-slope-at-inflection sigmoidmodels were fit to triplicate wells at seven concentrations to determinethe response characteristics for BIPI inhibitor compounds A, B and C.Each compound was evaluated using the full population of cellularmeasurements found in a 10 × 10-field scan area within each well, andalso limiting the number of cells measured to 100 in each well. R_(min)could not be established independently in compounds A, B (see FIGS.21a/21b and 23a/23b), and so R_(min) was fixed at 0.212 for thesecompounds (all compounds with fixed R_(min) marked with asterisks), asdetermined by the full population analysis of compound C. The G₁subpopulation response was also measured for compound B. Cells wereprepared as described above.

Analysis of Inhibitor Compounds: In another wellplate, cells weretreated with different concentrations of three inhibitor compounds (A, Band C), and then analyzed to assess the inhibition of NFκBtranslocation. In these experiments, both the full number of cells foundin a 10×10 scan area, and a fixed 100 cell set were measured andcompared. All three compounds responded clearly (FIGS. 21 a-23 b), withnearly identical uninhibited translocation (R_(max)) in the range0.3772-0.4239. Compounds A and B failed to plateau at maximum inhibitorconcentration, and response estimates failed to converge with a floatingR_(min) parameter. Compound C achieved maximum inhibition withR_(min)=0.212 (full cell population) and 0.220 (100 cell population);R_(min) was subsequently fixed to achieve response convergence oncompounds A and B. Table 13 summarizes all estimated response parameterscalculated by nonlinear optimization (unity slope at inflection isassumed in the response model). IC₅₀s varied 1.02×10⁻⁶-1.15×10⁻⁶ forcompound C, and by a factor in the approximate range 32-80 greater forcompounds A and B.

TABLE 14 Scan rate timing model components with definitions andestimated default values. tuning estimated component time required tovalue t_(plate, M×N) scan a M × N plate see equation t_(well, Q×R) scana Q × R field area inside a well see equation t_(field) image a singlefield see equation t_(I, k) integrate the k^(th) image channel ≧0.0167 st_(MW) move between adjacent wells and setup 1.00 s scan t_(MF) movebetween adjacent fields 0.30 s t_(RW) return from end of field rowinside a well 0.50 s t_(F) focus on a single field 0.40 s t_(T) transferimage data from camera 0.11 s t_(E) reconfigure emission/excitationfilters 0.20 s t_(S) open/close shutter 0.03 s t_(C) adjust cameraparameters 0.10 s

TABLE 15 System timing examples showing typical system throughputdependence on assay and cell preparation protocol, assuming two channelscans, 3 × 3 well scan areas, 96 well plates and estimated componenttimes listed in Table 14. t_(1, 2) t_(field) t_(well, 3×3)t_(plate, 8×12) 0.0167 s 0.9834 s 11.65 s 20.22 m 0.5000 s 1.4667 s 16.0 s 27.18 m A bright fluorescent intensity (t_(I, 2) = 0.0167 s)requires 0.9834 s/well and 20.22 m/plate, while a dimmer fluorescentintensity (t_(I, 2) = 0.5000 s), perhaps due to fluorophore binding to aless concentrated target, requires 16.0 s/well and 27.18 m/plate. Bothsecondary channels assume a bright nuclear fluorophore as the primarychannel (t_(I, 1) = 0.0167 s). Sufficient online image processing poweris assumed available so that throughput is scan-hardware limited.

Scan Rate Estimates: A simple model of scan rates is broken into plate,well and field components

$\begin{matrix}{t_{{plate},{M \times N}} = {{MNt}_{{well},{Q \times R}} + {\left( {{MN} - 1} \right)t_{M\; W}}}} & (19) \\{t_{{well},{Q \times R}} = {{QRt}_{field} + {{Q\left( {R - 1} \right)}t_{MF}} + {\left( {Q - 1} \right)t_{RW}}}} & (20) \\{t_{field} = {t_{S} + t_{F} + {\sum\limits_{k = 1}^{n}t_{I,k}} + {n\; t_{T}} + {\left( {n - 1} \right)\left( {t_{E} + t_{C}} \right)}}} & (21)\end{matrix}$

with definitions and estimated values given in Table 14. The assumptionsin this model include: (1) wellplates are scanned in a zigzag pattern,(2) wells are scanned in a raster pattern requiring a return betweenrows to maximize field alignment, (3) online processing occurscompletely in the background so that system is scan-hardware limited.

Table 15 gives timing example estimates for scanning a 2-channel, 96wellplate with 3×3 well scan areas for the two cases of a bright and dim(requiring integration) secondary channel (examples both assume a brightprimary nuclear channel), resulting in 11.65 s/well, 20.22 m/plate forthe bright assay, and 16.0 s/well, 27.18 m/plate for the dim assay.These examples are not specific to a particular application, which mayrequire more or less time depending on parameters, but are intended tosuggest the scope of scan rates that will be typical for assaysdeveloped for the platform of FIG. 1.

When comparing scan rates for the platform of FIG. 1 with other imagingplatforms, t_(F) and Q×R are the dominant factors. We found thatt_(F)=0.4 s was achieved due to the speed of the autofocus technologydescribed above. The well scan area Q×R is chosen to allow imaging of asample size large enough to accommodate some minimum significantresponse. Since a required minimum significant response (MSR) fixes n,we can define the ratio of the required field acquisition times of twosystems or assays as the time efficiency (TE) for comparison

$\begin{matrix}{{{TE}({MSR})} = \frac{{n_{A}({MSR})}d_{B}t_{{field},A}}{{n_{B}({MSR})}d_{A}t_{{field},B}}} & (22)\end{matrix}$

where A, B are comparable systems and d is the cellular field density(cells/field). TE accounts for differences in number of cells necessarybetween systems, but for simplicity, neglects times related to stagemotion.

Comparing the systems from Table 6 for NFκB translocation (B is theplatform of FIG. 1) at α=0.05, β=0.20, d_(B)/d_(A)=0.25, MSR=20% andletting t_(field,A)/t_(field,B) vary between 2-6 gives a TE range of5.95-17.85. Although this efficiency comparison relates strictly to NFκBtranslocation assays as reported here and in Ding et al., the rangeillustrates fundamental advantages over the platform described in Dinget al. afforded by the autofocus and image segmentation speed, and inimage measurement fidelity that will manifest in other assays (specificTEs would have to be calculated for each). The high end of the TE rangehere would become possible, for example, if the protocol were optimizedfor the platform of FIG. 1 to brighten the secondary channel, therebytaking better advantage of the autofocus speed. Scan rate comparisonsbetween NFκB assays can also be cast into an absolute time scale byadditionally assuming specific cells/field densities. FIG. 25 shows time(seconds) versus MSR (percent of total translocation dynamic range)curves for both the platform of FIG. 1 and the platform used by Ding etal., under typical experimental conditions.

FIG. 25 is a family of curves produced by plotting time required toachieve a minimum significant response ((MSR) as a percent of total NFκBtranslocation on maximum stimulation) per well. In these curves (lowestto highest), σ=0.040, best case for the platform of FIG. 1; σ=0.093,worst case for the platform of FIG. 1; σ=0.277=2.98×0.093, the platformof Ding et al., 2.98=ratio of measurement CVs (see Table 7).Assumptions: α=β=0.05, t_(field,Dingetal)=3 s,t_(field,figure1platform)=1.4667 s, d_(Dingetal)=40 cells/field (10×objective magnification), d_(figure1platform)=10 cells/field (20×objective magnification).

An Automatic Tool for Achieving Minimum Significant Response in a DrugScreen The error measurements described above can be used to predict thenumber of cells needed to reach a required minimum significant responsefor in a given compound screen. A graphical user interface tool willlead the user through the steps of measuring control and test wells todetermine the errors. This predictor tool will allow the user tomanipulate the time required to scan a plate as a function of theminimum significant response for the cell population and/orsubpopulation. With this automation, the use can easily determinewhether or not to use a subpopulation or the full population and howmuch time will be required to reach a given significance. The scan willthen proceed automatically by scanning each well until the requirednumber of cells is measured.

Refer now to FIGS. 26 a, 26 b, 27, 28, and 29 for an understanding of amethod of performing an automated drug screen according to theinvention. These figures include screens of a graphical user interface(GUI) designed to afford a user the ability to select and set parametersto perform the automated drug screen using the principles and functionsof this invention; they also include flow charts representinginitialization and performance of a drug screen according to theinvention. The illustrated method represents functionality invested inthe high throughput platform illustrated in FIG. 1 by programming thehost computer 110 with a series of software instructions. It may also besaid that these figures illustrate a method performable by a programmedcomputer. The GUI screens represent visible output produced on themonitor 115; they also represent functionality for identifying to a userparameters of a drug screen, and for enabling the user to input valuesfor those parameters to the host computer 110. These screens may also bereferred to as “GUI tools”.

Referring now to these figures, FIGS. 26 a, 26 b, and 28 showuser-executed steps carried out prior to performing the automated drugscreen in order to initialize parameters that control the screen. Theuser first chooses an assay that may be run for many days or weeks. NFκBis an example of such an assay. The Type I and Type II errors (equations8-16) are then entered into the host computer 110 through a GUI toolshown in FIG. 26 a. Control experiments are then performed to measurethe variability and amplitude of FLIC. When these values are at hand,they are entered into the host computer 110 through a GUI tool shown inFIG. 26 b. The variability is the standard deviation (SD) or coefficientof variation (CV) of FLIC and the percent amplitude is computed as100×(maximum-minimum)/(instrument measurement range). The maximum andminimum are determined by the cellular responses to controls (usuallywith no compound and with a maximally stimulating compound).

The automated drug screen is then performed by the high-throughputplatform of FIG. 1. Refer to FIGS. 27 and 29 for an understanding of howtis screen is performed. First, screening parameters are entered using aGUI tool such as that shown in FIG. 27. If limiting the time of the drugscreen is the primary concern, the user enters a desired time per plate(typically in minutes) using the Time/Plate field of the GUI tool, andthe estimated MSR and Cells/Well are calculated automatically anddisplayed. If the resulting MSR is deemed inappropriate, the user maychoose to shorten or lengthen the Time/Plate as required to achieve thedesired MSR. The MSR also depends on the number of cells per image whichis also entered by the user using the Cells/Image field of the GUI tool.The actual MSR will be a function of the number of cells per imageacquired as the drug screen is performed. If achieving a certain MSR isdeemed more important, the user enters the desired MSR using the MSRfield of the GUI tool and enters an average number of cells using theCells/Image field. Then the Cells/Well value and estimated Time/Plateare calculated automatically and displayed in the corresponding fieldsof the GUI tool. In this case, the actual time will depend on the valuein the Cells/Image field and the high throughput platform of FIG. 1 willacquire images until the desired MSR is achieved. The user can alsoenter the value in the Cells/Well field and the average number ofCells/Image and the MSR and Time/Plate will be calculated automaticallyand displayed in the corresponding fields.

The flow chart of FIG. 29 shows the steps executed by the highthroughput platform of FIG. 1 in automated scanning. The drug screen mayconsist of scanning a single microtiter plate (with 96, 384 or anynumber of wells) or a stack of microtiter plates delivered by a platerobot for a multi-plate scan. The scan is initiated in step 2900 by theuser after placing a plate on the microscope stage (or the plate isloaded by the robot). The platform then moves the plate to the firstwell (step 2910), performs autofocus (step 2920), acquires the colorimage in step 2930 (each color represents a different fluorescent dyestaining a different cellular molecule and/or cellular compartment),performs image segmentation in step 2940 (finds the masks of thecellular compartments), separates compartments of overlapping cellsusing tessellation, also in step 2940, and measures FLIC in step 2950.In step 2960 a test is then performed to determine if enough cells havebeen measured to achieve the desired MSR or if the number of images hasbeen acquired to achieve the desired time/plate. If not, in step 2970the platform moves to the next field of view (image) within the samewell and repeats the acquisition and analysis steps. When enough imageshave been analyzed, the platform moves the plate to the next well instep 2910 and repeats the same steps. When the last well in theexperiment has been analyzed (step 2980), the plate is completed andeither the next plate is automatically loaded (step 2990), or the scanis complete. When all of the plates in the experiment have beenanalyzed, the drug screen is finished.

System Design Considerations for Tessellation

An O(n log n) Voronoi Tessellation Algorithm: With reference to FIGS. 30a-30 h, this section describes the algorithm used to create Voronoitessellations for a collection of random points in the plane. The pointscome from an image analysis program that computes the position ofvarious objects from a scanning cytometer. The objects may be cellnuclei or other images. The only important fact is that the objects canbe represented by a set of nodes specified by X and Y coordinates. Forcells, cell nuclei and other objects the assignee's U.S. Pat. Nos.5,548,661 and 5,790,692, describe image segmentation techniques that canbe used to locate the objects in the image. Once those or other imagesegmentation techniques have been used to create a mask of pixelsdescribing the locations of the object pixels of each object, thecentroids or other single-pixel node can be computed that represents theobject for tessellation.

A useful objective is to tessellate a magnified, segmented image so thateach node lies inside its own polygon. The interior of the polygon iscloser to its node than to any other node. This tessellation is calledthe Voronoi tessellation, and it can be generated directly, or from itsdual graph, a Delaunay triangulation of the plane, in which the nodesform the vertices of triangles. Refer to FIG. 30 a, which shows a set ofnodes 3000 in a plane 3001. In the tessellation algorithm employed inour invention, the four corners 3002 of the plane are added asadditional nodes and the resulting set of nodes is triangulated, asshown in FIG. 30 b.

Many different triangulations are possible. The triangulation shown inFIG. 30 c has the “Delaunay property,” namely that the circle thatcircumscribes each triangle contains no nodes in its interior. In FIG.30 c, such a circle 3008 (called a “circumdisk”) circumscribing thecentral triangle contains no other nodes. Each of the triangles has itsown circumdisk, and no node lies in the interior of a circumdick

Note in FIG. 30 c that each node is the vertex of several triangles,which can be ordered either clockwise or counterclockwise. By connectingthe centers of the circumdisks of these triangles (for example, by line3010), one obtains the “Voronoi tessellation,” as shown in FIG. 30 d. AVoronoi tessellation is dual to the Delaunay triangulation. The polygonedges are perpendicular bisectors of the edges of the triangles. Theregion is now split up into polygons around each node. Each polygoncontains all points closer to its node than to any other node. Note thatthe four corner nodes also have polygons assigned to them.

Ours is an efficient algorithm to create the Delaunay triangulation andVoronoi tessellation for an arbitrary collection of nodes in the plane.According to our algorithm, if n is the number of nodes, then ouralgorithm scales linearly with n in the memory requirements, and scaleswith n log n for time requirements. This is close to optimal. Using ouralgorithm, we have found that 100000 nodes can be triangulated andpolygonized in under 5 seconds, with memory requirements on the order of35 MB. A million nodes can be done in around 75 seconds, with 379 MBrequired for storage. The machine used is a Wintel box running at 800MHz.

Presume computation of the Delaunay triangulation for a set of nodes.The minimal set of data needed for this is the following:

-   -   1) Node coordinates: An X and Y coordinate for each node, stored        as floating point numbers. The most convenient way to do this is        to assign each node an integer index, starting with 0 and going        up to n−1, where n is the number of nodes. Then the X and Y        coordinates are stored in two large double-precision arrays.    -   2) Triangle nodes: The indices of the three nodes that make up        each triangle. The easiest way to store this data is by        assigning each triangle an integer index, starting with 0 and        going up to N−1, where N is the number of triangles. N is        approximately 2n. Then one stores the nodes of each triangle in        three integer arrays, in counterclockwise order going around the        triangle.

It is also convenient (though not strictly necessary) to store thecenter coordinates and radius of the circumdisk of each triangle. Thesewill be stored in double-precision arrays indexed by the triangle index.Given the coordinates of each of the three nodes of a triangle, it is asimple linear algebra problem to compute the coordinates of the centerof the circumdisk. Then the radius is given by the Euclidean distancefrom the center to any one of the nodes.

We next describe an algorithm for efficiently finding the coordinatesfor all the polygons in the diagram in one pass. The algorithm is asfollows:

-   -   1) Create an array of lists of triangle indices, one list for        each node. So the array of lists will be indexed by the node        index.    -   2) For each triangle, get its three node indices. There is a        list for each of these nodes. To each of these three lists, add        the triangle index.    -   3) To construct the Voronoi polygon for any node, get its list        of triangles. Order the list counterclockwise by computing the        angle that the center of the circumdisk of each triangle makes        with the node, and then sorting in increasing order. The centers        of the circumdisks, ordered in this way, form the vertices of        the Voronoi polygon for the node.

All of this can be implemented easily using Standard Template Library(STL) components. The coordinates and triangle nodes can be stored invectors of doubles and integers. The array of lists can be implementedas a vector of vectors of integers. To minimize storage, the vectorsshould be sized in advance using the reserve( )method of the vectorclass. For example, the coordinate arrays should have n entriesreserved. The nodal arrays for triangles should have N entries reserved.The array of lists should have n lists reserved, with each list having 8entries reserved. It is rare for a node to have more than 8 triangles;in the few cases where a node has more than 8, the STL vector willautomatically resin itself. This algorithm scales in time linearly withthe number of nodes, and it executes much faster than the triangulationalgorithm.

We now describe an iterative method for triangulating a set of nodescontained in some rectangle. At each stage, a valid Delaunaytriangulation of k nodes will be assumed, and then the next node will beadded and edges will be rearranged so as to yield a Delaunaytriangulation of k+1 nodes.

Initially, a Delaunay triangulation of the four corner nodes of thecontaining triangle is constructed. This is shown in FIG. 30 e. With aDelaunay triangulation of the four corners of a rectangle, the samecircle circumscribes both triangles, and its center is at the midpointof the diagonal. Now presume successful triangulation a collection of knodes, and that we want to add the next node. Assume that the new nodelies in the interior of the existing triangulation. Since thedescription started with a bounding box, this is a safe assumption.

The algorithm for modifying the mesh is illustrated in FIG. 30 f:

-   -   1) Find the set of all triangles whose circumdisks contain the        new node. These triangles will have to be deleted, since a valid        Delaunay triangulation never has a triangle whose circumdisks        contain nodes.    -   2) The union of these triangles forms a polygon containing the        new node. Delete all the triangle edges contained in the        interior of this polygon.    -   3) Create new edges connecting the new node to each of the        vertices of this polygon.

The node 3050 is the new node being inserted. It is contained in thecircumdisks of four previously existing triangles, having a total of sixnodes 3051. The six exterior edges connecting the nodes 3051 are leftalone. The three interior edges shown by dashed lines are deleted. Sixnew edges (the black solid lines) are created, each connecting the newnode 3050 to one of the outer nodes. By construction the circumdisk ofeach of the new triangles intersects the new node, and the new node isnot in the interior of any other circumdisk. Hence, the newtriangulation is still Delaunay. The operation has increased the totalnode count by one, the total triangle count by two, and the total edgecount by three. This accounting holds true for each node added to thetriangulation. Therefore, the final mesh will have approximately 2ntriangles and 3n edges, where n is the number of nodes. Note that theinitial mesh has 4 nodes, 2 triangles, and 5 edges. Accordingly, theexact result is that the mesh will have 2n−6 triangles and 3n−7 edges.

Thus far, we have described the algorithm above only generally. Acomplete description requires a specification of the data structures andthe steps needed to carry out the steps. The simplest procedure is touse no new data structures. Step (1) can then be performed by abrute-force search of all triangles. Each triangle whose circumdiskcontains the new node is added to a list. In practice, this list issmall; it almost always contains fewer than 10 entries. Steps (2) and(3) can then be performed quickly. However, the brute-force search instep (1) requires an operation for each triangle, which is approximatelytwice the number of nodes. Since the brute-force approach requires asearch for each node, the time to triangulate the mesh requires of ordern² operations, and this is rather poor.

The algorithm thus far described is known as Watson's algorithm, and hasbeen around since the 1970s. In the next section, we describe animprovement to this algorithm that substantially speeds up thetriangulation process for large meshes.

The bottleneck in the triangulation scheme described above occurs whereall circumdisks containing the new node to be inserted are found. Ifeven one triangle whose circumdisk contains the new node can be found,then the search can be restricted to the nearest neighbors, which wouldsubstantially increase efficiency. So the problem reduces to findingthat first triangle whose circumdisk contains the new node.

Our preferred solution is use of the so-called “sweepline” method. Asweepline algorithm for computing the Voronoi tessellation was publishedby S. Fortune in 1987. Our algorithm uses something similar forcomputing the Delaunay triangulation. The idea is to sort the nodes inincreasing order of the X-coordinate. Ties would be broken by sorting onthe Y-coordinate. Since we assume that all (x, y) pairs are distinct,this defines a unique ordering of the nodes to be inserted. We theninsert nodes in this order, using the algorithm above, but now using anintelligent search rather than the brute-force search. One can imagine avertical line sweeping across from left to right. At any time, thevertical line is located at the X-coordinate most recently inserted. Allnodes to the left of the line have been inserted; all nodes to the rightare not inserted yet.

A brute force sweepline algorithm searches backwards in the list oftriangles. This list will have the most recently inserted triangles atthe tail of the list. If the algorithm starts at the tail of the list,it will search the right-most triangles and will likely find one quicklywhose circumdisk contains the next node to be inserted. This simplesolution works well, and reduces the triangulation time to order n logn. But there is a faster search algorithm, which is also reasonablysimple, so we'll describe that here.

A faster search algorithm depends on the fact that nodes are insertedinside a rectangle that is triangulated to begin with. So the right sideof the mesh consists of a large number of long, skinny triangles thatcontain one of the two corner nodes on the right, either the northeastor the southeast node. Our sweepline algorithm is illustrated in FIG. 30g. In this figure, the four corner nodes define the maximum extent ofthe region, and were triangulated first. The nodes 3062 in the interiorof the rectangle have been inserted from left to right. The shadedregion 3065 represents a large number (possibly millions) of nodes whichhave already been inserted. The vertical sweepline 3070 marks the lineof demarcation. To its left, all nodes have been inserted into the mesh.To its right, no nodes have been inserted. The red node is the next nodein the list. Since it lies to the right of the sweepline, it must lie inone of the long skinny triangles containing either the northeast node3060 n or the southeast node 3060 s (or it can lie in the large “centraltriangle” that contains both the northeast and southeast nodes). Thereare roughly sqrt(n) long skinny triangles connecting to the cornernodes. We put these into two lists, the “northeast list” and the“southeast list” and sort them by the slope of their lower edges. Thenthe node to be inserted will either be in the central triangle or it canbe rapidly located in one of these lists using a binary search. The timeto perform this search is of order log(sqrt(n)). It should be clear thatthe time to triangulate all the nodes is of order n log(sqrt(n)).Technically, this is still of order n log n, but it should be clear thatthe operation count is pretty small.

We have implemented this algorithm and its performance is of order n logn all the way up, to half a million nodes, (which can be triangulated inless than 30 seconds). We see slightly more time required for a millionnodes, but the reason is not clear. It could be memory allocationproblems, or some other subtlety. We anticipate the need to do mesheswith 100000 nodes as quickly as possible. The current implementation cancreate the mesh in under 5 seconds. We expect that this can be reducedby 10 to 20% by tweaking the code.

We now describe how to complete the list of triangles whose circumdiskcontains the node to be inserted once the first triangle has been foundusing the algorithm given above. The solution to this requiresbookkeeping. The idea is that all such triangles are adjacent. So assoon as the first one is found, one need only search all the nearestneighbors in some orderly way, until all of the possibilities have beeneliminated. This is implemented as a recursive operation.

First, we create an STL set containing the indices of all triangleswhose circumdisk contains the new node. We initialize it with the firsttriangle found, and then call a method that tests the triangle andrecursively checks its three neighbors. The recursion stops if atriangle is already in the set, or if its circumdisk doesn't contain thenode. We have found experimentally that these sets tend to be small andthat it is actually faster to use an STL vector, to avoid the costs ofinsertion into a set. The effect is a speed boost of about 12%.

Next the neighboring triangles of each triangle must be found. Thesimplest way to do this is to maintain a map of each edge to thetriangle on its right, as shown in FIG. 30 h. In this figure, The threenodes of the shaded triangle 3080, N₁, N₂, and N₃ define three edgesgoing counterclockwise around the triangle. The ordered pair (N₁, N₂) ismapped to the index of the triangle 3082 on their right. Likewise, theordered pair (N₂,N₁) is mapped to the shaded triangle 3080.

The easiest way to do this mapping is to define an STL map that takes apair of integers (the node indices) to a third integer (the index of thetriangle to the right). However, for large meshes, the memoryrequirements are quite large. Furthermore, this is relatively slow. Oneshould remember that there are a lot of edges. For each node in themesh, there are approximately 6 directed edges.

An easier way to do the mapping is to make an STL vector of maps. Thevector is indexed by nodes. The idea is that each node in the meshpoints to several others (on average, 6 others). So each entry in thevector will be a map of an integer to an integer, mapping a node indexto a triangle index. This is substantially more efficient than the mapof integer pairs to integers.

However, even better results are obtainable. The map of an integer to aninteger is still inefficient, since the main operations are insertingand deleting entries. It is faster and more space efficient to replacethem by vectors of pairs of integers. The first entry in the pair is anode index; the second entry in the pair is the triangle index. So ourdata structure is a vector of vectors of pairs of integers. It isimperative to use the RESERVE method (a function in Microsoft® VisualC++®) to allocate a reasonable amount of space for the vectors.Otherwise, the STL will allocate several kB for each vector, and thememory requirements will be outlandish.

This is a common theme in working with the STL. Small sets (or smallmaps) are not as efficient as small vectors (or small vectors of pairs).Sets and maps are far more convenient to work with for the programmer,but they just can't match the efficiency of a vector when you're dealingwith a handful of items. Ordinarily, this doesn't matter. But when onetries to optimize the inner loop, it can make sense to replace the setwith a vector, as long as a small amount of space is reserved for thevector. The STL vector will grow automatically if necessary, but it'sbest to make this a rare event. In our case, nodes typically average 6neighboring nodes, with 8 neighbors being unusual, and more than 8 beingrare.

1. A system for measuring cell activity represented in an image obtainedwith an automated microscopy platform the image including cells treatedwith an agent that causes translocation of cellular material, the systemcomprising: means for segmenting cellular components in the image; meansfor separating segmented components of overlapping cells in the image;and means for determining translocation of cellular material betweencellular compartments in the segmented, separated image based upon aratio of a first intensity of the cellular material in a first cellularcompartment to a sum of at least the first intensity and a secondintensity of the cellular material in a second cellular compartment. 2.The system of claim 1, wherein the means for segmenting includes meansfor enhancing contrast between the first cellular compartments and thebackground of the image, and between the second cellular compartmentsand the background of the image.
 3. The system of claim 1, wherein themeans for separating includes means for tessellating the image.
 4. Thesystem of claim 1, wherein the means for determining includes: means formeasuring a fractionalized local intensity of the first cellularmaterial; means for measuring a fractionalized local intensity of thesecond cellular material; and, means for comparing the fractionalizedlocal intensities of the first and second cellular materials.
 6. Thesystem of claim 4, wherein: the means for segmenting includes means forenhancing contrast between the first cellular compartments and thebackground of the image, and between the second cellular compartmentsand the background of the image; and, the means for separating includesmeans for tessellating the image.
 7. The system of claim 1, wherein themeans for determining includes: means for measuring the first intensityof the cellular material in the first cellular compartment a and formeasuring the second intensity of the cellular material in the secondcellular compartment b; and, means for determining a fractionallocalized intensity F_(a) of the cellular material in the first cellularcompartment by a ratio of the first intensity to a sum of at least thefirst and second intensities, according to:$F_{a} = \frac{\sum\limits_{{({x,y})} \in a}{I\left( {x,y} \right)}}{\sum\limits_{{{({x,y})} \in a},b}{I\left( {x,y} \right)}}$where I(x, y), the emission intensity at a pixel location (x, y) in acellular compartment, is determined according to:I(x, y) = ∫₀^(z)I₀Q ɛ uz in which I₀ is an incident excitationintensity, Q is a quantum yield, ε is an extinction coefficient, u is alocal fluorophore concentration at pixel location (x, y), and z is thethickness of a column of cellular material at pixel location (x, y). 8.The system of claim 7, wherein segmenting includes enhancing contrastbetween first cellular compartments including the first cellularcompartment and the background of the image, and between second cellularcompartments including the second cellular compartment and thebackground of the image.
 9. The system of claim 8, wherein separatingincludes tessellating the image.
 10. The system of claim 7, whereinseparating includes tessellating the image.